-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathanalytical.py
executable file
·122 lines (89 loc) · 3.71 KB
/
analytical.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
from numpy import *
from scipy import special
from scipy.special import lpmv
from scipy.misc import factorial
def get_K(x,n):
K = 0.
n_fact = factorial(n)
n_fact2 = factorial(2*n)
for s in range(n+1):
K += 2**s*n_fact*factorial(2*n-s)/(factorial(s)*n_fact2*factorial(n-s)) * x**s
return K
def solution(q, xq, E_1, E_2, R, kappa, a, N):
qe = 1.60217646e-19
Na = 6.0221415e23
E_0 = 8.854187818e-12
cal2J = 4.184
PHI = zeros(len(q))
for K in range(len(q)):
rho = sqrt(sum(xq[K]**2))
zenit = arccos(xq[K,2]/rho)
azim = arctan2(xq[K,1],xq[K,0])
phi = 0.+0.*1j
for n in range(N):
for m in range(-n,n+1):
P1 = lpmv(abs(m),n,cos(zenit))
Enm = 0.
for k in range(len(q)):
rho_k = sqrt(sum(xq[k]**2))
zenit_k = arccos(xq[k,2]/rho_k)
azim_k = arctan2(xq[k,1],xq[k,0])
P2 = lpmv(abs(m),n,cos(zenit_k))
Enm += q[k]*rho_k**n*factorial(n-abs(m))/factorial(n+abs(m))*P2*exp(-1j*m*azim_k)
C2 = (kappa*a)**2*get_K(kappa*a,n-1)/(get_K(kappa*a,n+1) +
n*(E_2-E_1)/((n+1)*E_2+n*E_1)*(R/a)**(2*n+1)*(kappa*a)**2*get_K(kappa*a,n-1)/((2*n-1)*(2*n+1)))
C1 = Enm/(E_2*E_0*a**(2*n+1)) * (2*n+1)/(2*n-1) * (E_2/((n+1)*E_2+n*E_1))**2
if n==0 and m==0:
Bnm = Enm/(E_0*R)*(1/E_2-1/E_1) - Enm*kappa*a/(E_0*E_2*a*(1+kappa*a))
else:
Bnm = 1./(E_1*E_0*R**(2*n+1)) * (E_1-E_2)*(n+1)/(E_1*n+E_2*(n+1)) * Enm - C1*C2
phi += Bnm*rho**n*P1*exp(1j*m*azim)
PHI[K] = real(phi)/(4*pi)
C0 = qe**2*Na*1e-3*1e10/(cal2J)
E_P = 0.5*C0*sum(q*PHI)
return E_P
def solution_2(q, xq, E_1, E_2, R, kappa, a, N):
qe = 1.60217646e-19
Na = 6.0221415e23
E_0 = 8.854187818e-12
cal2J = 4.184
PHI = zeros(len(q))
for K in range(len(q)):
rho = sqrt(sum(xq[K]**2))
zenit = arccos(xq[K,2]/rho)
azim = arctan2(xq[K,1],xq[K,0])
phi = 0.+0.*1j
for n in range(N):
for m in range(-n,n+1):
P1 = lpmv(abs(m),n,cos(zenit))
Enm = 0.
for k in range(len(q)):
rho_k = sqrt(sum(xq[k]**2))
zenit_k = arccos(xq[k,2]/rho_k)
azim_k = arctan2(xq[k,1],xq[k,0])
P2 = lpmv(abs(m),n,cos(zenit_k))
# Enm += q[k]*rho_k**n*factorial(n-abs(m))/factorial(n+abs(m))*P2*exp(-1j*m*azim_k)
Enm += q[k]*rho_k**n*conjugate(special.sph_harm(m, n, azim_k, zenit_k))
C2 = (kappa*a)**2*get_K(kappa*a,n-1)/(get_K(kappa*a,n+1) +
n*(E_2-E_1)/((n+1)*E_2+n*E_1)*(R/a)**(2*n+1)*(kappa*a)**2*get_K(kappa*a,n-1)/((2*n-1)*(2*n+1)))
C1 = Enm/(E_2*E_0*a**(2*n+1)) * (2*n+1)/(2*n-1) * (E_2/((n+1)*E_2+n*E_1))**2
if n==0 and m==0:
Bnm = (Enm/(E_0*R)*(1/E_2-1/E_1) - Enm*kappa*a/(E_0*E_2*a*(1+kappa*a)))*4*pi/(2*n+1)
else:
Bnm = (1./(E_1*E_0*R**(2*n+1)) * (E_1-E_2)*(n+1)/(E_1*n+E_2*(n+1)) * Enm - C1*C2)* 4*pi/(2*n+1)
# phi += Bnm*rho**n*P1*exp(1j*m*azim)
phi += Bnm*rho**n*special.sph_harm(m, n, azim, zenit)
PHI[K] = real(phi)/(4*pi)
C0 = qe**2*Na*1e-3*1e10/(cal2J)
E_P = 0.5*C0*sum(q*PHI)
return E_P
#q = array([1.])
#xq = array([[1e-12,1e-12,1e-12]])
#E_1 = 4.
#E_2 = 80.
#R = 1.
#kappa = 0.125
#a = 1.
#N = 20
#print(solution(q, xq, E_1, E_2, R, kappa, a, N))
#print(solution_2(q, xq, E_1, E_2, R, kappa, a, N))