在 Python 中从 IDL 定义函数

Defining a function from IDL in Python

我正在尝试将大量代码从 IDL 复制到 Python。我认为我遇到的一个问题是函数 gridgen 的定义。 Gridgen 是一个函数,用于生成具有相等对数间距网格点的垂直网格,其中: zmin = 网格顶部的坐标; zmax = 网格底部的坐标; Nlvls = 网格中所需的级别数,z = 输出网格。

IDL代码是:

FUNCTION GRIDGEN, zmin, zmax, Nlvls
  dlnz       = (ALOG(zmax)-ALOG(zmin))/(Nlvls-1) ;divisions in log space
  z          = FLTARR(Nlvls) ;array of Nlvls points for logarithm to base 10
  z[*]       = 0.            ;initialize grid to zero
  z[Nlvls-1] = ALOG(zmax)    ;assign the maximum value in log spacing
  FOR i=Nlvls-2, 0, -1 DO z[i] = z[i+1] - dlnz ;generate log spacing values
  z          = EXP(z) ;convert from log(z) to actual values
  RETURN, z
END

我如何将其翻译成 Python 是:

def gridgen100(zmin, zmax, Nlvls):
    dlnz = ((np.log(zmax) - np.log(zmin))/(Nlvls - 1))  # divisions in log space
    z = np.zeros(Nlvls, dtype=float)                    # array of Nlvls points for logarithm to base 10
    z[Nlvls-1] = np.log(zmax)                           # assign the maximum value in log spacing
    for i in np.arange(Nlvls-2, Nlvls-101, -1):         # NOT CORRECT; correct is: for i in [Nlvls-2, 0, -1]:
      z[i] = z[i +1] - dlnz                             # generate log spacing values
      #z = np.exp(np.array(z))                          # convert from log(z) to actual values [MUST DO OUTSIDE DEF]
return z

问题是:

  1. 由于我使用 np.arange 创建 for 循环的方式,当我有不同的 Nlvls 时,我必须定义一个单独的 gridgen 函数(例如,有时我有 100 个 Nlvls,有时我有 25 个)。
  2. 我无法在函数内将 log(z) 转换为实际值,我必须在定义之外进行。

我目前无法访问 IDL,因此我无法通过比较 IDL 输出与 Python 输出来排除故障。

我是 Python 自学成才的初学者,但我很感激任何人可以提供的任何帮助或建议。

IIUC,您的 IDL for 循环转换为

for i in range(Nlvls-2, -1, -1):

即从 Nlvls-2 开始,递减 1 直到达到 0。这给了我

def gridgen(zmin, zmax, Nlvls):
    dlnz = (np.log(zmax) - np.log(zmin))/(Nlvls-1)
    z = np.zeros(Nlvls, dtype=float)
    z[Nlvls-1] = np.log(zmax)
    for i in range(Nlvls-2, -1, -1):
        z[i] = z[i+1] - dlnz
    z = np.exp(z)
    return z

>>> gridgen(2, 8, 10)
array([ 2.        ,  2.33305808,  2.72158   ,  3.1748021 ,  3.70349885,
        4.32023896,  5.0396842 ,  5.87893797,  6.85795186,  8.        ])

但是已经有一个 numpy 函数,np.logspace,它会为你做这个日志间距,所以如果我对你所追求的是正确的,你可以使用它来获得相同的结果:

>>> np.logspace(np.log10(2), np.log10(8), 10)
array([ 2.        ,  2.33305808,  2.72158   ,  3.1748021 ,  3.70349885,
        4.32023896,  5.0396842 ,  5.87893797,  6.85795186,  8.        ])

(对于第 2 点,如果不想 return 到原来的 space,显然可以删除 z = np.exp(z) 行。)