Rebecca M. Rebecca M. - 7 months ago 11
Python Question

Defining a function from IDL in Python

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.

DSM DSM
Answer

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