Python函数返回类型与fortran子例程声明类型不匹配

2024-09-28 05:17:31 发布

您现在位置:Python中文网/ 问答频道 /正文

我有一些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:

^{pr2}$

并在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)

任何帮助都将不胜感激。在


Tags: selfk2k1nnrdngdoenddo

热门问题