1 """
2 Calculate the stellar surface displacements due to pulsations.
3
4
5 """
6 import numpy as np
7 from numpy import sqrt,pi,sin,cos,exp
8 from scipy.special import lpmv,legendre
9 from scipy.misc.common import factorial
10 from scipy.integrate import dblquad,quad
11 from scipy.spatial import Delaunay
12
13
14
16 """
17 Legendre polynomial.
18
19 Check equation (3) from Townsend, 2002:
20
21 >>> ls,x = [0,1,2,3,4,5],cos(linspace(0,pi,100))
22 >>> check = 0
23 >>> for l in ls:
24 ... for m in range(-l,l+1,1):
25 ... Ppos = legendre_(l,m,x)
26 ... Pneg = legendre_(l,-m,x)
27 ... mycheck = Pneg,(-1)**m * factorial(l-m)/factorial(l+m) * Ppos
28 ... check += sum(abs(mycheck[0]-mycheck[1])>1e-10)
29 >>> print check
30 0
31 """
32 m_ = abs(m)
33 legendre_poly = legendre(l)
34 deriv_legpoly_ = legendre_poly.deriv(m=m_)
35 deriv_legpoly = np.polyval(deriv_legpoly_,x)
36 P_l_m = (-1)**m_ * (1-x**2)**(m_/2.) * deriv_legpoly
37 if m<0:
38 P_l_m = (-1)**m_ * factorial(l-m_)/factorial(l+m_) * P_l_m
39 return P_l_m
40
42 """
43 Spherical harmonic according to Townsend, 2002.
44
45 This function is memoized: once a spherical harmonic is computed, the
46 result is stored in memory
47
48 >>> theta,phi = mgrid[0:pi:20j,0:2*pi:40j]
49 >>> Ylm20 = sph_harm(theta,phi,2,0)
50 >>> Ylm21 = sph_harm(theta,phi,2,1)
51 >>> Ylm22 = sph_harm(theta,phi,2,2)
52 >>> Ylm2_2 = sph_harm(theta,phi,2,-2)
53 >>> p = figure()
54 >>> p = gcf().canvas.set_window_title('test of function <sph_harm>')
55 >>> p = subplot(411);p = title('l=2,m=0');p = imshow(Ylm20.real,cmap=cm.RdBu)
56 >>> p = subplot(412);p = title('l=2,m=1');p = imshow(Ylm21.real,cmap=cm.RdBu)
57 >>> p = subplot(413);p = title('l=2,m=2');p = imshow(Ylm22.real,cmap=cm.RdBu)
58 >>> p = subplot(414);p = title('l=2,m=-2');p = imshow(Ylm2_2.real,cmap=cm.RdBu)
59
60 """
61 factor = (-1)**m * sqrt( (2*l+1)/(4*pi) * factorial(l-m)/factorial(l+m))
62 Plm = legendre_(l,m,cos(theta))
63 return factor*Plm*exp(1j*m*phi)
64
66 """
67 Derivative of spherical harmonic wrt colatitude.
68
69 Using Y_l^m(theta,phi).
70
71 Equation::
72
73 sin(theta)*dY/dtheta = (l*J_{l+1}^m * Y_{l+1}^m - (l+1)*J_l^m * Y_{l-1,m})
74
75 E.g.: Phd thesis of Joris De Ridder
76 """
77 if abs(m)>=l:
78 Y = 0.
79 else:
80 factor = 1./sin(theta)
81 term1 = l * norm_J(l+1,m) * sph_harm(theta,phi,l+1,m)
82 term2 = (l+1) * norm_J(l,m) * sph_harm(theta,phi,l-1,m)
83 Y = factor * (term1 - term2)
84 return Y/sin(theta)
85
87 """
88 Derivative of spherical harmonic wrt longitude.
89
90 Using Y_l^m(theta,phi).
91
92 Equation::
93
94 dY/dphi = i*m*Y
95 """
96 return 1j*m*sph_harm(theta,phi,l,m)
97
98
100 """
101 Normalisation factor
102 """
103 if abs(m)<l:
104 J = sqrt( float((l**2-m**2))/(4*l**2-1))
105 else:
106 J = 0
107 return J
108
110 """
111 Omega is actually spin parameter (Omega_rot/omega_freq)
112 """
113 return Omega * (l-abs(m)+1)/(l+1) * 2./(2*l+1) * (1-l*k)
114
116 """
117 Omega is actually spin parameter (Omega_rot/omega_freq)
118 """
119 return Omega * (l+abs(m))/l * 2./(2*l+1) * (1 + (l+1)*k)
120
121
122
123
124
126 """
127 Radial displacement, see Zima 2006.
128
129 t in phase units
130 """
131 return sph_harm(theta,phi,l,m) * exp(1j*t)
132
134 term1 = k * dsph_harm_dtheta(theta,phi,l,m) * exp(1j*t)
135 term2 = norm_atlp1(l,m,Omega,k) / sin(theta) * dsph_harm_dphi(theta,phi,l+1,m) * exp(1j*t + pi/2)
136 term3 = norm_atlm1(l,m,Omega,k) / sin(theta) * dsph_harm_dphi(theta,phi,l-1,m) * exp(1j*t - pi/2)
137 return term1 + term2 + term3
138
140 term1 = k /sin(theta) * dsph_harm_dphi(theta,phi,l,m)*exp(1j*t)
141 term2 = -norm_atlp1(l,m,Omega,k) * dsph_harm_dtheta(theta,phi,l+1,m)*exp(1j*t+pi/2)
142 term3 = -norm_atlm1(l,m,Omega,k) * dsph_harm_dtheta(theta,phi,l-1,m)*exp(1j*t-pi/2)
143 return term1 + term2 + term3
144
145 -def surface(theta,phi,l,m,t,Omega=0.1,k=1.,asl=0.2,radius=1.):
146 ksi_r = asl*sqrt(4*pi)*radial(theta,phi,l,m,t)
147 if l>0:
148 ksi_theta = asl*sqrt(4*pi)*colatitudinal(theta,phi,l,m,t,Omega,k)
149 ksi_phi = asl*sqrt(4*pi)*longitudinal(theta,phi,l,m,t,Omega,k)
150 else:
151 ksi_theta = np.zeros_like(theta)
152 ksi_phi = np.zeros_like(phi)
153
154 return (radius+ksi_r.real),\
155 (theta + ksi_theta.real),\
156 (phi + ksi_phi.real)
157
158 -def observables(theta,phi,teff,logg,l,m,t,Omega=0.1,k=1.,asl=0.2,radius=1.,delta_T=0.01+0j,delta_g=0.01+0.5j):
159 rad_part = radial(theta,phi,l,m,t)
160 ksi_r = asl*sqrt(4*pi)*rad_part
161 if l>0:
162 ksi_theta = asl*sqrt(4*pi)*colatitudinal(theta,phi,l,m,t,Omega,k)
163 ksi_phi = asl*sqrt(4*pi)*longitudinal(theta,phi,l,m,t,Omega,k)
164 else:
165 ksi_theta = np.zeros_like(theta)
166 ksi_phi = np.zeros_like(phi)
167 gravity = 10**(logg-2)
168 return (radius+ksi_r.real),\
169 (theta + ksi_theta.real),\
170 (phi + ksi_phi.real),\
171 (teff + (delta_T*rad_part*teff).real),\
172 np.log10(gravity+(delta_g*rad_part*gravity).real)+2
173
174
175
176 if __name__=="__main__":
177 from ivs.roche import local
178 from ivs.coordinates import vectors
179 from enthought.mayavi import mlab
180 from divers import multimedia
181 theta,phi = local.get_grid(50,25,gtype='triangular')
182 r = np.ones_like(theta)
183 x,y,z = vectors.spher2cart_coord(r,phi,theta)
184 points = np.array([x,y,z]).T
185 grid = Delaunay(points)
186
187
188 l,m = 2,2
189 asl = 0.01
190
191 for k in [0,1.,2.]:
192 for l in range(1,5):
193 for m in range(0,l+1,1):
194 mlab.figure(size=(1000,800))
195 mlab.gcf().scene.disable_render = True
196 if l==0 or l==1:
197 asl = 0.1
198 else:
199 asl = 0.01
200 old_center=None
201 for i,t in enumerate(np.linspace(0,2*pi,100)):
202 print k,l,m,i
203 r,th,ph = surface(theta,phi,l,m,t,asl=asl,k=k)
204 center,size,normal = local.surface_normals(r,ph,th,grid,gtype='triangular')
205
206 if i==0:
207 colors = r
208 r_c,phi_c,theta_c = vectors.cart2spher_coord(*center.T)
209 colors_ = r_c
210
211
212 mlab.clf()
213 mlab.points3d(center.T[0],center.T[1],center.T[2],colors_,scale_factor=0.05,scale_mode='none',colormap='RdBu',vmin=colors_.min(),vmax=colors_.max())
214
215 mlab.colorbar()
216
217 if i>=1:
218 vx,vy,vz = center.T[0]-old_center.T[0],\
219 center.T[1]-old_center.T[1],\
220 center.T[2]-old_center.T[2]
221 v = np.sqrt(vx**2+vy**2+vz**2)
222 mlab.quiver3d(center.T[0],center.T[1],center.T[2],\
223 vx,vy,vz,scalars=v,colormap='spectral',scale_mode='scalar')
224
225 old_center = center.copy()
226 mlab.view(distance=5,azimuth=-90,elevation=90)
227 mlab.savefig('pulsation_lm%d%d_k%03d_%03d.png'%(l,m,k,i))
228 mlab.close()
229 multimedia.make_movie('pulsation_lm%d%d_k%03d_*.png'%(l,m,k),output='pulsation_lm%d%d_k%03d.avi'%(l,m,k))
230