code: def get_Sb(data):
      
      Sb, N = 0., 0.
      for rec in data:
          delta = 90. - abs(rec['vgp_lat'])
          if rec['average_k'] != 0:
              k = rec['average_k']
              L = rec['average_lat'] * np.pi / 180.  # latitude in radians
              Nsi = rec['average_nn']
              K = old_div(k, (2. * (1. + 3. * np.sin(L)**2) /
                              (5. - 3. * np.sin(L)**2)))
              Sw = old_div(81., np.sqrt(K))
          else:
              Sw, Nsi = 0, 1.
          Sb += delta**2. - old_div((Sw**2), Nsi)
          N += 1.
      return np.sqrt(old_div(Sb, float(N - 1.)))
docs: returns vgp scatter for data set