code: def PG_ISSO_eq(r, chi, ci):
      
      X=chi**2*(chi**2*(3*chi**2+4*r*(2*r-3))+r**2*(15*r*(r-4)+28))-6*r**4*(r**2-4)
      Y=chi**4*(chi**4+r**2*(7*r*(3*r-4)+36))+6*r*(r-2)*(chi**6+2*r**3*(chi**2*(3*r+2)+3*r**2*(r-2)))
      Z=ISCO_eq(r, chi)
  
      return r**8*Z+chi**2*(1-ci**2)*(chi**2*(1-ci**2)*Y-2*r**4*X)
docs: Polynomial that enables the calculation of a generic innermost
      stable spherical orbit (ISSO) radius via its roots.  Physical
      solutions are between the equatorial ISSO (aka the ISCO) radius
      and the polar ISSO radius.
      [See Stone, Loeb, Berger, PRD 87, 084053 (2013)]
  
      Parameters
      -----------
      r: float
          the radial coordinate in BH mass units
      chi: float
          the BH dimensionless spin parameter
      ci: float
          cosine of the inclination angle between the BH spin
          and the orbital angular momentum
  
      Returns
      ----------
      float
          r**8*Z+chi**2*(1-ci**2)*(chi**2*(1-ci**2)*Y-2*r**4*X)
          where
          X=chi**2*(chi**2*(3*chi**2+4*r*(2*r-3))+r**2*(15*r*(r-4)+28))-6*r**4*(r**2-4)
          Y=chi**4*(chi**4+r**2*(7*r*(3*r-4)+36))+6*r*(r-2)*(chi**6+2*r**3*(chi**2*(3*r+2)+3*r**2*(r-2)))
          Z=ISCO_eq(r,chi)