我有一些Python代码如下,spherical.py
:
...
import f_utils
# (r, theta, phi)
class spherical_utilities(object):
def __init__(self, ng, n, lab):
self.n = n
self.ng = ng
self.set_ngauss()
self.set_funcs_ang()
self.set_all_layers(lab)
def set_ngauss(self):
k, w = ps_roots(self.ng)
self.knots, self.weights = k * pi, w * pi
self.thetas = self.knots
self.sint = sin(self.thetas)
self.cost = cos(self.thetas)
self.tgt = tan(self.thetas)
self.ctgt = 1 / self.tgt
def set_funcs_ang(self):
P = array([special.lpmn(self.n, self.n, ct) for ct in self.cost])
self.data_Ang, self.data_Angd = P[:, 0, :, :], P[:, 1, :, :]
def set_all_layers(self, lab):
self.data_layers = [0]
for bnd in lab.boundaries():
r, rd, rdd = bnd.shape.R(self.knots)
rdr = rd / r
r2rd2 = r ** 2 + rd ** 2
Rad = [0, {}, {}]
Radd = [0, {}, {}]
kis = [0, bnd.k1, bnd.k2]
for i in [1, 2]:
krs = kis[i] * r
JY = array([special.sph_jnyn(self.n, kr) for kr in krs])[:, :, :]
Rad[i] = {'j': JY[:, 0, :], 'h': JY[:, 0, :] + 1j * JY[:, 2, :]}
Radd[i] = {'j': JY[:, 1, :], 'h': JY[:, 1, :] + 1j * JY[:, 3, :]}
self.data_layers.append({'ki': kis,\
'r': r,\
'rd': rd,\
'rdd': rdd,\
'rdr': rdr,\
'r2rd2': r2rd2,\
'Rad': Rad,\
'Radd': Radd})
...
def mat_ebcm_axi_tm(C, m, jh1, jh2, e12, e21):
k1 = C.ki[1]
k2 = C.ki[2]
Rad1 = C.Rad(m, jh1, 1)
Radd1 = C.Radd(m, jh1, 1)
Rad2 = C.Rad(m, jh2, 2)
Radd2 = C.Radd(m, jh2, 2)
Angm = C.Ang(m)
Angmd = C.Angd(m)
return f_utils.axitm(Rad1, Radd1, Rad2, Radd2, Angm, Angmd,\
C.r, C.rd, C.sint, C.cost, k1, k2, e12, C.weights)
在其他文件中,methods.py
:
并在f_utils.for
中有一些对应的fortran代码:
SUBROUTINE axitm (NG,NN,Rad1,Radd1,Rad2,Radd2,Ang,Angd,
& Rs,Rsd,SinT,CosT,k1,k2,ep12,w,O)
implicit none
INTEGER NG,NN
COMPLEX*16 Rad1(NG,NN),Radd1(NG,NN),Rad2(NG,NN),Radd2(NG,NN)
COMPLEX*16 Ang(NG,NN), Angd(NG,NN)
COMPLEX*16 k1,k2,ep12
REAL*8 SinT(NG),CosT(NG), Rs(NG),Rsd(NG)
COMPLEX*16 O(NN,NN)
COMPLEX*16 w(NG)
cf2py intent(out) O
cf2py intent(hide) NG,NN
integer l,n,k
complex*16 coef, j1,jd1,j2,jd2,pl,pdl,pn,pdn
real*8 sn,cs,r,rd
complex*16 fun1
DO l=1,NN
DO n=1,NN
O(l,n)=0d0
ENDDO
ENDDO
DO k=1,NG
r = Rs(k)
rd = Rsd(k)
cs = CosT(k)
sn = SinT(k)
DO l=1,NN
DO n=1,NN
j1 = Rad1(k,l)
jd1 = Radd1(k,l)
j2 = Rad2(k,n)
jd2 = Radd2(k,n)
pl = Ang(k,l)
pdl = Angd(k,l)
pn = Ang(k,n)
pdn = Angd(k,n)
fun1 = k1**2 * r**2 * ( jd1 * j2 -
& ep12 * k2/k1 * j1 * jd2) *
& pl * pn * sn
& + k1 * rd * sn**2 * (pdl * pn - ep12*
& pl * pdn) * j1 * j2
& - (ep12 - 1) * (k1 * r * sn - k1 * rd * cs) *
& j1 * j2 * pl * pn
O(l,n) = O(l,n) + fun1*w(k)
ENDDO
ENDDO
ENDDO
DO l=1,NN
coef = (0d0,1d0) * (2d0 * l + 1) / (2 * l * (l+1))
DO n=1,NN
O(l,n) = O(l,n) * coef
ENDDO
ENDDO
end
对于Python:
return f_utils.axitm(Rad1, Radd1, Rad2, Radd2, Angm, Angmd,\
C.r, C.rd, C.sint, C.cost, k1, k2, e12, C.weights)
对于fortran:
SUBROUTINE axitm (NG,NN,Rad1,Radd1,Rad2,Radd2,Ang,Angd,
& Rs,Rsd,SinT,CosT,k1,k2,ep12,w,O)
我不想编辑fortran和Python的原始文件。在
当我正常运行spherical.py
时,我得到以下错误:
ValueError: 0-th dimension must be fixed to 1 but got 200
对于Python中的这一行代码:
return f_utils.axitm(Rad1, Radd1, Rad2, Radd2, Angm, Angmd,\
C.r, C.rd, C.sint, C.cost, k1, k2, e12, C.weights)
任何帮助都将不胜感激。在
目前没有回答
相关问题 更多 >
编程相关推荐