1 """
2 Roche models of fast and differential rotation
3
4 Section 1. Non-rotating spherically symmetric single star
5 =========================================================
6
7 The most trivial usage of this module would be to construct a non-rotating,
8 spherically symmetric single star. All quantities should be constant over the
9 stellar surface, except of course the projected ones. We compute this via the
10 zero-rotation limit of the fast rotating model:
11
12 First set some parameters: the rotation rate, effective temperature at the pole,
13 polar radius, mass and viewing angle.
14
15 >>> omega = 0.00 # ratio to critical velocity
16 >>> T_pole = 5500. # K
17 >>> r_pole = 1. # solar radii
18 >>> M = 1. # solar mass
19 >>> view_angle = pi/2 # radians
20
21 Construct a coordinate grid with a resolution of 50 points, covering the whole
22 star. Unravel, so that we have access to 1D coordinate arrays
23
24 >>> theta,phi = local.get_grid(50,full=True)
25 >>> thetas,phis = np.ravel(theta),np.ravel(phi)
26
27 Calculate the shape of the stellar surface and the local gravity for every point
28 in the grid. This is a 3D vector in Cartesian coordinates, so reshape all
29 components to match the grid shape. As a reference, also explicitly calculate
30 the polar surface gravity, which is the z-component of the gravity vector.
31
32 >>> radius = np.array([get_fastrot_roche_radius(itheta,r_pole,omega) for itheta in thetas]).reshape(theta.shape)
33 >>> 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
34 >>> grav_local = np.array([i.reshape(theta.shape) for i in grav_local])
35 >>> g_pole = fastrot_roche_surface_gravity(r_pole,0,0,r_pole,omega,M)[-1]
36
37 Calculate the value of the surface gravity in SI-units:
38
39 >>> grav = vectors.norm(grav_local)
40
41 Compute the size of all surface elements, and the angle with respect to the
42 radius vector (should be all zero). The normal to the surface is the negative of
43 the surface gravity vector, which is directed inwards:
44
45 >>> areas_local,cos_gamma = local.surface_elements((radius,(theta,phi)),-grav_local)
46
47 Now we can calculate the local effective temperature (assuming radiative
48 atmosphere in Von Zeipel law), and bolometric intensities:
49
50 >>> teff_local = local.temperature(vectors.norm(grav_local),g_pole,T_pole,beta=1.)
51 >>> ints_local = local.intensity(teff_local,grav,photband='OPEN.BOL')
52
53 Define the line-of-sight vector:
54
55 >>> line_of_sight = np.zeros_like(grav_local)
56 >>> line_of_sight[0] = -sin(view_angle)
57 >>> line_of_sight[2] = -cos(view_angle)
58
59 Compute the angles between surface normals and line-of-sight. We do this in 1D:
60
61 >>> gravity = np.array([igrav.ravel() for igrav in grav_local])
62 >>> line_of_sight = np.array([ilos.ravel() for ilos in line_of_sight])
63 >>> angles = vectors.angle(-gravity,line_of_sight)
64 >>> mus = cos(angles)
65
66 Now compute the bolometric projected intensity, taking care of limbdarkening
67 effects:
68
69 >>> intens_proj = local.intensity(teff_local.ravel(),grav.ravel(),mus,photband='OPEN.BOL').reshape(theta.shape)
70
71 The total and projected luminosity can then be calculated the following way:
72
73 >>> print pi*(ints_local*areas_local*constants.Rsol_cgs**2).sum()/constants.Lsol_cgs
74 0.992471247895
75 >>> print pi*np.nansum(intens_proj*areas_local*constants.Rsol_cgs**2)/constants.Lsol_cgs
76 0.360380373413
77
78 Now make some plots showing the local quantities:
79
80 >>> quantities = areas_local,np.log10(grav*100),teff_local,ints_local,intens_proj,angles/pi*180,radius
81 >>> names = 'Area','log g', 'Teff', 'Flux', 'Proj. flux', 'Angle'
82 >>> p = pl.figure()
83 >>> rows,cols = 2,3
84 >>> for i,(quantity,name) in enumerate(zip(quantities,names)):
85 ... p = pl.subplot(rows,cols,i+1)
86 ... p = pl.title(name)
87 ... q = quantity.ravel()
88 ... vmin,vmax = q[-np.isnan(q)].min(),q[-np.isnan(q)].max()
89 ... 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()
90 ... p = pl.scatter(phis/pi*180,thetas/pi*180,c=q,edgecolors='none',vmin=vmin,vmax=vmax)
91 ... p = pl.colorbar()
92 ... p = pl.xlim(0,360)
93 ... p = pl.ylim(180,0)
94
95 ]include figure]]ivs_binary_rochepotential_nonrotstar.png]
96
97 Section 2. Fast and uniformly rotating star
98 ===========================================
99
100 Repeating the same example as above, but now with parameters
101
102 >>> omega = 0.98
103
104 Gives the figure below:
105
106 ]include figure]]ivs_binary_rochepotential_fastrotstar.png]
107
108 Changing the viewing angle to 80 degrees
109
110 >>> view_angle = 80/180.*pi
111
112 gives: Note that the maximum projected intensity is higher than in the previous
113 example: this is because we view directly at higher latitudes now, which are
114 hotter than the equator. The least flux is coming from the equatorial zone.
115
116 ]include figure]]ivs_binary_rochepotential_fastrotstar_80incl.png]
117
118 To generate an image of the star, it is wise to convert the coordinates to
119 Cartesian coordinates:
120
121 >>> x,y,z = vectors.spher2cart_coord(radius.ravel(),phis,thetas)
122 >>> p = pl.figure()
123 >>> p = pl.subplot(121,aspect='equal')
124 >>> p = pl.title('top view')
125 >>> p = pl.scatter(x,y,c=teff_local.ravel(),edgecolors='none')
126 >>> p = pl.subplot(122,aspect='equal')
127 >>> p = pl.title('side view')
128 >>> p = pl.scatter(y,z,c=teff_local.ravel(),edgecolors='none')
129
130 The movie below is obtained by calculating the projected intensity in the line
131 of sight
132
133 >>> proj_intens,mus = local.projected_intensity(teff_local,grav_local,areas_local,np.array([0,-sin(view_angle),-cos(view_angle)]),photband='OPEN.BOL')
134
135 for different lines of sight, and then for each line of sight computing the
136 Cartesian coordinates, and finally rotating in Y-Z plane.
137
138 >>> y,z = vectors.rotate(y,z,-(pi/2-view_angle))
139
140 ]include figure]]ivs_binary_rochepotential_fastrotstar_shape.png]
141
142 ]include figure]]ivs_binary_rochepotential_fastrot.gif]
143
144 Section 3. Differentially rotating star
145 =======================================
146
147 The differential rotation rate has to follow the stereotypical law
148
149 Omega(theta) = Omega_pole + b * sin(theta)
150
151 The critical velocity can now easily be exceeded at the equator, if the rotation
152 rate increases towards the pole.
153
154 First set some parameters: the rotation rate of the equator and pole, effective
155 temperature at the pole, polar radius, mass and viewing angle.
156
157 >>> omega_eq = 0.1 # ratio to critical rotation rate (equatorial)
158 >>> omega_pl = 0.9 # ratio to critical rotation rate (polar)
159 >>> T_pole = 5500. # K
160 >>> r_pole = 1.0 # solar radii
161 >>> M = 1. # solar mass
162 >>> view_angle = pi/2 # radians
163
164 We now do very similar stuff as in Section 1, except for the different Roche
165 potential. (We can skip making the grid now)
166
167 >>> radius = np.array([get_diffrot_roche_radius(itheta,r_pole,M,omega_eq,omega_pl) for itheta in thetas]).reshape(theta.shape)
168 >>> 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
169 >>> grav_local = np.array([i.reshape(theta.shape) for i in grav_local])
170 >>> g_pole = diffrot_roche_surface_gravity(r_pole,0,0,r_pole,M,omega_eq,omega_pl)[-1]
171
172 Calculate the value of the surface gravity in SI-units:
173
174 >>> grav = vectors.norm(grav_local)
175
176 Compute the size of all surface elements, and the angle with respect to the
177 radius vector (should be all zero). The normal to the surface is the negative of
178 the surface gravity vector, which is directed inwards:
179
180 >>> areas_local,cos_gamma = local.surface_elements((radius,(theta,phi)),-grav_local)
181
182 Now we can calculate the local effective temperature (assuming radiative
183 atmosphere in Von Zeipel law), and bolometric intensities:
184
185 >>> teff_local = local.temperature(vectors.norm(grav_local),g_pole,T_pole,beta=1.)
186 >>> ints_local = local.intensity(teff_local,grav,photband='OPEN.BOL')
187
188 Compute the angles between surface normals and line-of-sight. We do this in 1D:
189
190 >>> gravity = np.array([igrav.ravel() for igrav in grav_local])
191 >>> line_of_sight = np.array([ilos.ravel() for ilos in line_of_sight])
192 >>> angles = vectors.angle(-gravity,line_of_sight)
193 >>> mus = cos(angles)
194
195 Now compute the bolometric projected intensity, taking care of limbdarkening
196 effects:
197
198 >>> intens_proj = local.intensity(teff_local.ravel(),grav.ravel(),mus,photband='OPEN.BOL').reshape(theta.shape)
199
200 Now make some plots showing the local quantities:
201
202 >>> quantities = areas_local,np.log10(grav*100),teff_local,ints_local,intens_proj,angles/pi*180,radius
203 >>> names = 'Area','log g', 'Teff', 'Flux', 'Proj. flux', 'Angle'
204 >>> p = pl.figure()
205 >>> rows,cols = 2,3
206 >>> for i,(quantity,name) in enumerate(zip(quantities,names)):
207 ... p = pl.subplot(rows,cols,i+1)
208 ... p = pl.title(name)
209 ... q = quantity.ravel()
210 ... vmin,vmax = q[-np.isnan(q)].min(),q[-np.isnan(q)].max()
211 ... 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()
212 ... p = pl.scatter(phis/pi*180,thetas/pi*180,c=q,edgecolors='none',vmin=vmin,vmax=vmax)
213 ... p = pl.colorbar()
214 ... p = pl.xlim(0,360)
215 ... p = pl.ylim(180,0)
216
217 ]include figure]]ivs_binary_rochepotential_diffrotstar.png]
218
219 To generate an image of the star, it is wise to convert the coordinates to
220 Cartesian coordinates:
221
222 >>> x,y,z = vectors.spher2cart_coord(radius.ravel(),phis,thetas)
223 >>> p = pl.figure()
224 >>> p = pl.subplot(121,aspect='equal')
225 >>> p = pl.title('top view')
226 >>> p = pl.scatter(x,y,c=teff_local.ravel(),edgecolors='none')
227 >>> p = pl.subplot(122,aspect='equal')
228 >>> p = pl.title('side view')
229 >>> p = pl.scatter(y,z,c=teff_local.ravel(),edgecolors='none')
230
231 ]include figure]]ivs_binary_rochepotential_diffrotstar_shape.png]
232 """
233 import numpy as np
234 from numpy import pi,cos,sin,sqrt,nan
235 from scipy.optimize import newton
236
237 from ivs.roche import local
238 from ivs.coordinates import vectors
239 from ivs.units import constants
240 from ivs.units import conversions
241
242
243
244
246 """
247 Calculate components of the local surface gravity of the fast rotating
248 Roche model.
249
250 Input units are solar units.
251 Output units are SI.
252 Omega is fraction of critical velocity.
253
254 See Cranmer & Owocki, Apj (1995)
255
256 @param r: radius of the surface element to calculate the surface gravity
257 @type r: float/ndarray
258 @param theta: colatitude of surface element
259 @type theta: float/ndarray
260 @param phi: longitude of surface element
261 @type phi: float/ndarray
262 @param r_pole: polar radius of model
263 @type r_pole: float
264 @param omega: fraction of critical rotation velocity
265 @type omega: float
266 @param M: mass of the star
267 @type M: float
268 @param norm: compute magnitude of surface gravity (True) or vector (False)
269 @type norm: boolean
270 @return: surface gravity magnitude or vector
271 @rtype: 3Xfloat/3Xndarray
272 """
273 GG = constants.GG_sol
274
275 x = r/r_pole
276 grav_r = GG*M/r_pole**2 * (-1./x**2 + 8./27.*x*omega**2*sin(theta)**2)
277
278 grav_th = GG*M/r_pole**2 * (8./27.*x*omega**2*sin(theta)*cos(theta))
279
280 grav = np.array([grav_r,grav_th])
281
282 grav = np.array(vectors.spher2cart( (r,phi,theta),(grav[0],0.,grav[1]) ))
283 grav = grav*constants.Rsol
284 if norm:
285 return vectors.norm(grav)
286 else:
287 return grav
288
290 """
291 Calculate Roche radius for a fast rotating star.
292
293 @param theta: angle from rotation axis
294 @type theta: float
295 @param r_pole: polar radius in solar units
296 @type r_pole: float
297 @param omega: angular velocity (in units of the critical angular velocity)
298 @omega_: float
299 @return: radius at angle theta in solar units
300 @rtype: float
301 """
302
303 Rstar = 3*r_pole/(omega*sin(theta)) * cos((pi + np.arccos(omega*sin(theta)))/3.)
304
305 if np.isinf(Rstar) or sin(theta)<1e-10:
306 Rstar = r_pole
307 return Rstar
308
310 """
311 Compute the critical angular velocity (Hz).
312
313 Definition taken from Cranmer and Owocki, 1995 and equal to
314
315 Omega_crit = sqrt( 8GM / 27Rp**3 )
316
317 Example usage (includes conversion to period in days):
318
319 >>> Omega = critical_angular_velocity(1.,1.)
320 >>> P = 2*pi/Omega
321 >>> P = conversions.convert('s','d',P)
322 >>> print 'Critical rotation period of the Sun: %.3f days'%(P)
323 Critical rotation period of the Sun: 0.213 days
324
325 @param M: mass (solar masses)
326 @type M: float
327 @param R_pole: polar radius (solar radii)
328 @type R_pole: float
329 @param units: if you wish other units than Hz, you can give them here
330 @type units: string understandable by L{<units.conversions.convert>}
331 @return: critical angular velocity in Hz
332 @rtype: float
333 """
334 M = M*constants.Msol
335 R = R_pole*constants.Rsol
336 omega_crit = np.sqrt( 8*constants.GG*M / (27.*R**3))
337 if units.lower()!='hz':
338 omega_crit = conversions.convert('Hz',units,omega_crit)
339 return omega_crit
340
342 """
343 Compute the critical velocity (km/s)
344
345 Definition 1 from Cranmer and Owocki, 1995:
346
347 v_c = 2 pi R_eq(omega_c) * omega_c
348
349 Definition 2 from Townsend 2004:
350
351 v_c = sqrt ( 2GM/3Rp )
352
353 which both amount to the same value:
354
355 >>> critical_velocity(1.,1.,definition=1)
356 356.71131858379499
357 >>> critical_velocity(1.,1.,definition=2)
358 356.71131858379488
359
360 @param M: mass (solar masses)
361 @type M: float
362 @param R_pole: polar radius (solar radii)
363 @type R_pole: float
364 @param units: if you wish other units than Hz, you can give them here
365 @type units: string understandable by L{units.conversions.convert}
366 @return: critical velocity in km/s
367 @rtype: float
368 """
369 if definition==1:
370 omega_crit = critical_angular_velocity(M,R_pole)
371 P = 2*pi/omega_crit
372 R_eq = get_fastrot_roche_radius(pi/2.,R_pole,1.)*constants.Rsol
373 veq = 2*pi*R_eq/P
374 elif definition==2:
375 veq = np.sqrt( 2*constants.GG * M*constants.Msol / (3*R_pole*constants.Rsol))
376 veq = conversions.convert('m/s',units,veq)
377
378 return veq
379
380
381
382
383
385 """
386 Definition of Roche potential due to differentially rotating star
387
388 We first solve the cubic equation
389
390 M{re/rp = 1 + f (x^2 + x + 1)/(6x^2)}
391
392 where M{f = re^3 Omega_e^2 / (G M)}
393 and M{x = Omega_e / Omega_p}
394
395 This transforms to solving
396
397 M{re^3 + b * re + c = 0}
398
399 where M{b = -1 / (aXrp) and c = 1/(aX),}
400 and M{a = Omega_e^2/(GM) and X = (x^2 + x + 1)/(6x^2)}
401
402 @param r: radius of the surface element to calculate the surface gravity
403 @type r: float/ndarray
404 @param theta: colatitude of surface element
405 @type theta: float/ndarray
406 @param r_pole: polar radius of model
407 @type r_pole: float
408 @param M: mass of the star
409 @type M: float
410 @param omega_eq: fraction of critical rotation velocity at equator
411 @type omega_eq: float
412 @param omega_pole: fraction of critical rotation velocity at pole
413 @type omega_pole: float
414 @return: roche potential value
415 @rtype: float/ndarray
416 """
417 GG = constants.GG_sol
418
419 Omega_crit = sqrt(8*GG*M/ (27*r_pole**3))
420 omega_eq = omega_eq*Omega_crit
421 omega_pole = omega_pole*Omega_crit
422 x = omega_eq / omega_pole
423
424
425 a = omega_eq**2/(GG*M)
426 X = (x**2+x+1)/(6*x**2)
427 b,c = -1./(a*X*r_pole), +1./(a*X)
428 m,k = 27*c+0*1j,-3*b+0*1j
429 n = m**2-4*k**3 + 0*1j
430 om1 = -0.5 + 0.5*sqrt(3)*1j
431 om2 = -0.5 - 0.5*sqrt(3)*1j
432 c1 = (0.5*(m+sqrt(n)))**(1./3.)
433 c2 = (0.5*(m-sqrt(n)))**(1./3.)
434 x1 = -1./3. * ( c1 + c2 )
435 x2 = -1./3. * ( om2*c1 + om1*c2 )
436 x3 = -1./3. * ( om1*c1 + om2*c2 )
437 re = x2.real
438
439
440 f = re**3 * omega_eq**2 / (GG*M)
441
442 rat = 1 + (f*(x**2+x+1))/(6.*x**2)
443
444 alpha = f*(x-1)**2/(6*x**2)*(1/rat)**7
445 beta = f*(x-1) /(2*x**2)*(1/rat)**5
446 gamma = f /(2*x**2)*(1/rat)**3
447
448 sinth = sin(theta)
449 y = r/r_pole
450 surf = alpha*y**7*sinth**6 + beta*y**5*sinth**4 + gamma*y**3*sinth**2 - y +1
451 return surf
452
454 """
455 Surface gravity from differentially rotation Roche potential.
456
457 Magnitude is OK, please carefully check direction.
458
459 @param r: radius of the surface element to calculate the surface gravity
460 @type r: float/ndarray
461 @param theta: colatitude of surface element
462 @type theta: float/ndarray
463 @param phi: longitude of surface element
464 @type phi: float/ndarray
465 @param r_pole: polar radius of model
466 @type r_pole: float
467 @param M: mass of the star
468 @type M: float
469 @param omega_eq: fraction of critical rotation velocity at equator
470 @type omega_eq: float
471 @param omega_pole: fraction of critical rotation velocity at pole
472 @type omega_pole: float
473 @param norm: compute magnitude of surface gravity (True) or vector (False)
474 @type norm: boolean
475 @return: surface gravity magnitude or vector
476 @rtype: 3Xfloat/3Xndarray
477 """
478 GG = constants.GG_sol
479 Omega_crit = sqrt(8*GG*M/ (27*r_pole**3))
480 omega_eq = omega_eq*Omega_crit
481 omega_pole = omega_pole*Omega_crit
482 x = omega_eq / omega_pole
483
484
485 a = omega_eq**2/(GG*M)
486 X = (x**2+x+1)/(6*x**2)
487 b,c = -1./(a*X*r_pole), +1./(a*X)
488 m,k = 27*c+0*1j,-3*b+0*1j
489 n = m**2-4*k**3 + 0*1j
490 om1 = -0.5 + 0.5*sqrt(3)*1j
491 om2 = -0.5 - 0.5*sqrt(3)*1j
492 c1 = (0.5*(m+sqrt(n)))**(1./3.)
493 c2 = (0.5*(m-sqrt(n)))**(1./3.)
494 x1 = -1./3. * ( c1 + c2 )
495 x2 = -1./3. * ( om2*c1 + om1*c2 )
496 x3 = -1./3. * ( om1*c1 + om2*c2 )
497 re = x2.real
498
499
500 f = re**3 * omega_eq**2 / (GG*M)
501
502 rat = 1 + (f*(x**2+x+1))/(6.*x**2)
503
504 alpha = f*(x-1)**2/(6*x**2)*(1/rat)**7
505 beta = f*(x-1) /(2*x**2)*(1/rat)**5
506 gamma = f /(2*x**2)*(1/rat)**3
507
508 sinth = sin(theta)
509 y = r/r_pole
510 grav_th = (6*alpha*y**7*sinth**5 + 4*beta*y**5*sinth**3 + 2*gamma*y**3*sinth)*cos(theta)
511 grav_r = 7*alpha/r_pole*y**6*sinth**6 + 5*beta/r_pole*y**4*sinth**4 + 3*gamma/r_pole*y**2*sinth**2 - 1./r_pole
512
513 fr = 6*alpha*y**7*sinth**4 + 4*beta*y**5*sinth**2 + 2*gamma*y**3
514 magn = GG*M/r**2 * sqrt(cos(theta)**2 + (1-fr)**2 *sin(theta)**2)
515 magn_fake = np.sqrt(grav_r**2+(grav_th/r)**2)
516 grav_r,grav_th = grav_r/magn_fake*magn,(grav_th/r)/magn_fake*magn
517
518 grav = np.array([grav_r*constants.Rsol,grav_th*constants.Rsol])
519
520 grav = vectors.spher2cart( (r,phi,theta),(grav[0],0.,grav[1]) )
521
522 if norm:
523 return vectors.norm(grav)
524 else:
525 return grav
526
527
529 """
530 Calculate Roche radius for a differentially rotating star.
531
532 @param theta: angle from rotation axis
533 @type theta: float
534 @param r_pole: polar radius in solar units
535 @type r_pole: float
536 @param M: mass in solar units
537 @type M: float
538 @param omega_eq: equatorial angular velocity (in units of the critical angular velocity)
539 @omega_eq: float
540 @param omega_pole: polar angular velocity (in units of the critical angular velocity)
541 @omega_pole: float
542 @return: radius at angle theta in solar units
543 @rtype: float
544 """
545 try:
546 r = newton(diffrot_roche_potential,r_pole,args=(theta,r_pole,M,omega_eq,omega_pole))
547 except RuntimeError:
548 r = np.nan
549 return r
550
552 """
553 Evaluate a differential rotation law of the form Omega = b1+b2*omega**2
554
555 The relative differential rotation rate is the ratio of the rotational shear
556 to the equatorial velocity
557
558 alpha = omega_eq - omega_pole / omega_eq
559
560 The units of angular velocity you put in, you get out (i.e. in terms of
561 the critical angular velocity or not).
562
563 @param omega_eq: equatorial angular velocity
564 @type omega_eq: float
565 @param omega_pole: polar angular velocity
566 @type omega_pole: float
567 @param theta: colatitude (0 at the pole) in radians
568 @type theta: float (radians)
569 @return: the angular velocity at colatitude theta
570 @rtype: float
571 """
572 b1 = omega_eq
573 b2 = omega_pole - omega_eq
574 omega_theta = b1 + b2 * cos(theta)**2
575 return omega_theta
576
577
579 """
580 Calculate the velocity vector of every surface element.
581
582 @param coordinates: polar coordinates of stellar surface (phi,theta,radius)
583 make sure the radius is in SI units!
584 @type coordinates: 3xN array
585 @param omega_eq: equatorial angular velocity (as a fraction of the critical one)
586 @type omega_eq: float
587 @param omega_pole: polar angular velocity (as a fraction of the critical one)
588 @type omega_pole: float
589 @param R_pole: polar radius in solar radii
590 @type R_pole: float
591 @param M: stellar mass in solar mass
592 @type M: float
593 @return: velocity vectors of all surface elements
594 @rtype 3xN array
595 """
596
597 Omega_crit = critical_angular_velocity(M,R_pole)
598 phi,theta,radius = coordinates
599 omega_local = diffrot_law(omega_eq,omega_pole,theta)*Omega_crit
600
601 omega_local_vec = np.array([np.zeros_like(omega_local),np.zeros_like(omega_local),omega_local]).T
602
603 x,y,z = vectors.spher2cart_coord(radius,phi,theta)
604 surface_element = np.array([x,y,z]).T
605
606 velo_local = np.array([np.cross(ielement,iomega_local_vec) for ielement,iomega_local_vec in zip(surface_element,omega_local_vec)]).T
607 return velo_local
608
609
610
611
612 if __name__=="__main__":
613 import doctest
614 import pylab as pl
615 doctest.testmod()
616 pl.show()
617