code: def _phi2deriv(self,R,z,phi=0.,t=0.):
          
          if not self.isNonAxi:
              phi= 0.
          x,y,z= self._compute_xyz(R,phi,z,t)
          Fx= self._xforce_xyz(x,y,z)
          Fy= self._yforce_xyz(x,y,z)
          Fxy= np.dot(self.rot(t, transposed = True),np.array([Fx,Fy]))
          Fx, Fy= Fxy[0], Fxy[1]
          phixxa= self._2ndderiv_xyz(x,y,z,0,0)
          phixya= self._2ndderiv_xyz(x,y,z,0,1)
          phiyya= self._2ndderiv_xyz(x,y,z,1,1)
          ang = self._omegab*t + self._pa
          c, s = np.cos(ang), np.sin(ang)
          phixx = c**2*phixxa + 2.*c*s*phixya + s**2*phiyya
          phixy = (c**2-s**2)*phixya + c*s*(phiyya - phixxa)
          phiyy = s**2*phixxa - 2.*c*s*phixya + c**2*phiyya
          return R**2.*(np.sin(phi)**2.*phixx+np.cos(phi)**2.*phiyy\
                            -2.*np.cos(phi)*np.sin(phi)*phixy)\
                            +R*(np.cos(phi)*Fx+np.sin(phi)*Fy)
docs: NAME:
             _phi2deriv
          PURPOSE:
             evaluate the second azimuthal derivative for this potential
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             the second azimuthal derivative