Defining a function from IDL in Python

Question

I am trying to copy an extensive code from IDL into Python. One issue I believe I am having is with the definition of the function gridgen. Gridgen is a function used to generate a vertical grid with equal log-spaced grid points, where: zmin = coordinate at top of grid; zmax = coordinate at bottom of grid; Nlvls = desired number of levels in grid and z = output grid.

The IDL code is:

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

How I translated that into Python is:

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

The issues are:

1. Because of how I created the for loop using np.arange, I have to define a separate gridgen function when I have different Nlvls (for example, sometimes I have 100 Nlvls, sometimes I have 25).
2. I can't convert from log(z) to actual values within the function, I have to do it outside the definition.

I don't currently have access to IDL, so I'm unable to troubleshoot by comparing the IDL output to the Python output.

I am self-taught in Python and a beginner, but I appreciate any help or advice anyone can offer.

Answer 1

IIUC, your IDL for loop translates to

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

i.e. start at Nlvls-2 and drop by 1 until you reach 0. This gives me

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

and

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

But there's already a numpy function, np.logspace, which does this log spacing for you, so if I'm right about what you're after, you could obtain the same result using it instead:

>>> 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.        ])

(For point #2, you can obviously remove the z = np.exp(z) line if you don't want to return to the original space.)