code: def epade(np, ns, ip, k0, dr, pd1, pd2):
      m = 40
      dg = range(m)
      dh1 = range(m)
      dh3 = range(m)
      a = [[0 for i in range(m)] for j in range(m)]
      b = range(m)
      bin = [[0 for i in range(m)] for j in range(m)]
      fact = range(m)
      ci = complex(0.0, 1.0)
      sig = k0 * dr
      n = 2 * np
      if ip == 1:
          nu = 0.0
          alp = 0.0
      else:
          nu = 1.0
          alp = -0.25
      #     The factorials
      fact[0] = 1.0
      for i in range(1, n):
          fact[i] = (i + 1) * fact[i - 1]
      #     The binomial coefficients.;
      for i in range(0, n + 1):
          bin[0][i] = 1.0
          bin[i][i] = 1.0
      for i in range(2, n + 1):
          for j in range(1, i):
              bin[j][i] = bin[j - 1][i - 1] + bin[j][i - 1]
      for i in range(0, n):
          for j in range(0, n):
              a[i][j] = 0.0
      #     The accuracy constraints.;
      deriv(n, sig, alp, dg, dh1, dh3, bin, nu)
      for i in range(0, n):
          b[i] = dg[i + 1]
      for i in range(1, n + 1):
          if 2 * i - 1 <= n:
              a[2 * i - 2][i - 1] = fact[i - 1]
          for j in range(1, i + 1):
              if 2 * j <= n:
                  a[2 * j - 1][i - 1] = (
                      -bin[j][i] * fact[j - 1] * dg[i - j]
                  )
      #     The stability constraints.;
      if ns >= 1:
          z1 = -3.0
          b[n - 1] = -1.0
          for j in range(1, np + 1):
              a[2 * j - 2][n - 1] = z1**j
              a[2 * j - 1][n - 1] = 0.0
      if ns >= 2:
          z1 = -1.5
          b[n - 2] = -1.0
          for j in range(1, np + 1):
              a[2 * j - 2][n - 2] = z1**j
              a[2 * j - 1][n - 2] = 0.0
      gauss(n, a, b)
      dh1[0] = 1.0
      for j in range(1, np + 1):
          dh1[j] = b[2 * j - 2]
      dh2 = fndrt(dh1, np)
      for j in range(0, np):
          pd1[j] = -1.0 / dh2[j]
      dh1[0] = 1.0
      for j in range(1, np + 1):
          dh1[j] = b[2 * j - 1]
      dh2 = fndrt(dh1, np)
      for j in range(0, np):
          pd2[j] = -1.0 / dh2[j]
  
file_path: pythran/pythran/tests/cases/ramsurf.py
imports: import cmath
label_desc: Use image node_id 6 for calling a global function with example usage: epade(np, ns, ip, k0, dr, pd1, pd2) without return types