code: def pOparapar(self,Opar,apar,tdisrupt=None):
          
          if _APY_LOADED and isinstance(Opar,units.Quantity):
              Opar= Opar.to(1/units.Gyr).value\
                  /bovy_conversion.freq_in_Gyr(self._vo,self._ro)
          if _APY_LOADED and isinstance(apar,units.Quantity):
              apar= apar.to(units.rad).value
          if tdisrupt is None: tdisrupt= self._tdisrupt
          if isinstance(Opar,(int,float,numpy.float32,numpy.float64)):
              Opar= numpy.array([Opar])
          out= numpy.zeros(len(Opar))
          # Compute ts
          ts= apar/Opar
          # Evaluate
          out[(ts < tdisrupt)*(ts >= 0.)]=\
              numpy.exp(-0.5*(Opar[(ts < tdisrupt)*(ts >= 0.)]-self._meandO)**2.\
                             /self._sortedSigOEig[2])/\
                             numpy.sqrt(self._sortedSigOEig[2])
          return out
docs: NAME:
  
             pOparapar
  
          PURPOSE:
  
             return the probability of a given parallel (frequency,angle) offset pair
  
          INPUT:
  
             Opar - parallel frequency offset (array) (can be Quantity)
  
             apar - parallel angle offset along the stream (scalar) (can be Quantity)
  
          OUTPUT:
  
             p(Opar,apar)
  
          HISTORY:
  
             2015-12-07 - Written - Bovy (UofT)