code: def computeORFMatrix(psr):
      
  
      # begin loop over all pulsar pairs and calculate ORF
      npsr = len(psr)
      ORF = N.zeros((npsr, npsr))
      phati = N.zeros(3)
      phatj = N.zeros(3)
      ptheta = [N.pi/2 - p['DECJ'].val for p in psr]
      pphi = [p['RAJ'].val for p in psr]
      for ll in range(0, npsr):
          phati[0] = N.cos(pphi[ll]) * N.sin(ptheta[ll])
          phati[1] = N.sin(pphi[ll]) * N.sin(ptheta[ll])
          phati[2] = N.cos(ptheta[ll])
  
          for kk in range(0, npsr):
              phatj[0] = N.cos(pphi[kk]) * N.sin(ptheta[kk])
              phatj[1] = N.sin(pphi[kk]) * N.sin(ptheta[kk])
              phatj[2] = N.cos(ptheta[kk])
             
              if ll != kk:
                  xip = (1.-N.sum(phati*phatj)) / 2.
                  ORF[ll, kk] = 3.*( 1./3. + xip * ( N.log(xip) -1./6.) )
              else:
                  ORF[ll, kk] = 2.0
  
      return ORF
docs: Compute ORF matrix.
  
      :param psr: List of pulsar object instances
  
      :returns: Matrix that has the ORF values for every pulsar
               pair with 2 on the diagonals to account for the 
               pulsar term.