code: def Rz_to_lambdanu(R,z,ac=5.,Delta=1.):
      
      g = Delta**2 / (1.-ac**2)
      a = g - Delta**2
      term  =  R**2 + z**2 - a - g
      discr = (R**2 + z**2 - Delta**2)**2 + (4. * Delta**2 * R**2)
      l = 0.5 * (term + nu.sqrt(discr))  
      n = 0.5 * (term - nu.sqrt(discr))
      if isinstance(z,float) and z == 0.:
          l = R**2 - a
          n = -g
      elif isinstance(z,nu.ndarray) and nu.sum(z == 0.) > 0:
          if isinstance(R,float):      l[z==0.] = R**2 - a
          if isinstance(R,sc.ndarray): l[z==0.] = R[z==0.]**2 - a
          n[z==0.] = -g
      return (l,n)
docs: NAME:
  
         Rz_to_lambdanu
  
      PURPOSE:
  
         calculate the prolate spheroidal coordinates (lambda,nu) from
         galactocentric cylindrical coordinates (R,z)            
         by solving eq. (2.2) in Dejonghe & de Zeeuw (1988a) for (lambda,nu):
              R^2 = (l+a) * (n+a) / (a-g)
              z^2 = (l+g) * (n+g) / (g-a)
              Delta^2 = g-a
  
      INPUT:
  
          R     - Galactocentric cylindrical radius
          z     - vertical height
          ac    - axis ratio of the coordinate surfaces 
                  (a/c) = sqrt(-a) / sqrt(-g) (default: 5.)
          Delta - focal distance that defines the spheroidal coordinate system (default: 1.)
                  Delta=sqrt(g-a)
  
      OUTPUT:
  
         (lambda,nu)
  
      HISTORY:
  
         2015-02-13 - Written - Trick (MPIA)