code: def pOparapar(self,Opar,apar):
          
          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 isinstance(Opar,(int,float,numpy.float32,numpy.float64)):
              Opar= numpy.array([Opar])
          out= numpy.zeros(len(Opar))
          # Compute ts and where they were at impact for all
          ts= apar/Opar
          apar_impact= apar-Opar*self._timpact
          dOpar_impact= self._kick_interpdOpar(apar_impact)
          Opar_b4impact= Opar-dOpar_impact
          # Evaluate the smooth model in the two regimes:
          # stripped before or after impact
          afterIndx= (ts < self._timpact)*(ts >= 0.)
          out[afterIndx]=\
              super(streamgapdf,self).pOparapar(Opar[afterIndx],
                                                apar)
          out[True^afterIndx]=\
              super(streamgapdf,self).pOparapar(Opar_b4impact[True^afterIndx],
                                                apar_impact[True^afterIndx],
                                                tdisrupt=
                                                self._tdisrupt-self._timpact)
          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-11-17 - Written - Bovy (UofT)