docs: def plot_eq_cont(fignum, DIblock, color_map='coolwarm'):
      """
      plots dec inc block as a color contour
      Parameters
      __________________
      Input:
          fignum :  figure number
          DIblock : nested pairs of [Declination, Inclination]
          color_map : matplotlib color map [default is coolwarm]
      Output:
          figure
      """
      import random
      plt.figure(num=fignum)
      plt.axis("off")
      XY = []
      centres = []
      counter = 0
      for rec in DIblock:
          counter = counter + 1
          X = pmag.dir2cart([rec[0], rec[1], 1.])
          # from Collinson 1983
          R = old_div(np.sqrt(1. - X[2]), (np.sqrt(X[0]**2 + X[1]**2)))
          XY.append([X[0] * R, X[1] * R])
      # radius of the circle
      radius = (old_div(3., (np.sqrt(np.pi * (9. + float(counter)))))) + 0.01
      num = 2. * (old_div(1., radius))  # number of circles
      # a,b are the extent of the grids over which the circles are equispaced
      a1, a2 = (0. - (radius * num / 2.)), (0. + (radius * num / 2.))
      b1, b2 = (0. - (radius * num / 2.)), (0. + (radius * num / 2.))
      # this is to get an array (a list of list wont do) of x,y values
      xlist = np.linspace(a1, a2, int(np.ceil(num)))
      ylist = np.linspace(b1, b2, int(np.ceil(num)))
      X, Y = np.meshgrid(xlist, ylist)
      # to put z in the array I just multiply both x,y with zero.  I will add to
      # the zero values later
      Z = X * Y * 0.
      # keeping the centres of the circles as a separate list instead of in
      # array helps later
      for j in range(len(ylist)):
          for i in range(len(xlist)):
              centres.append([xlist[i], ylist[j]])
      # the following lines are to figure out what happens at the edges where part of a circle might lie outside
      # a thousand random numbers are generated within the x,y limit of the circles and tested whether it is contained in
      # the eq area net space....their ratio gives the fraction of circle
      # contained in the net
      fraction = []
      beta, alpha = 0.001, 0.001  # to avoid those 'division by float' thingy
      for i in range(0, int(np.ceil(num))**2):
          if np.sqrt(((centres[i][0])**2) + ((centres[i][1])**2)) - 1. < radius:
              for j in range(1, 1000):
                  rnd1 = random.uniform(
                      centres[i][0] - radius, centres[i][0] + radius)
                  rnd2 = random.uniform(
                      centres[i][1] - radius, centres[i][1] + radius)
                  if ((centres[i][0] - rnd1)**2 + (centres[i][1] - rnd2)**2) <= radius**2:
                      if (rnd1**2) + (rnd2**2) < 1.:
                          alpha = alpha + 1.
                          beta = beta + 1.
                      else:
                          alpha = alpha + 1.
              fraction.append(old_div(alpha, beta))
              alpha, beta = 0.001, 0.001
          else:
              fraction.append(1.)  # if the whole circle lies in the net
  
      # for every circle count the number of points lying in it
      count = 0
      dotspercircle = 0.
      for j in range(0, int(np.ceil(num))):
          for i in range(0, int(np.ceil(num))):
              for k in range(0, counter):
                  if (XY[k][0] - centres[count][0])**2 + (XY[k][1] - centres[count][1])**2 <= radius**2:
                      dotspercircle += 1.
              Z[i][j] = Z[i][j] + (dotspercircle * fraction[count])
              count += 1
              dotspercircle = 0.
      im = plt.imshow(Z, interpolation='bilinear', origin='lower',
                      # cmap=plt.color_map.hot, extent=(-1., 1., -1., 1.))
                      cmap=color_map, extent=(-1., 1., -1., 1.))
      plt.colorbar(shrink=0.5)
      x, y = [], []
      # Draws the border
      for i in range(0, 360):
          x.append(np.sin((old_div(np.pi, 180.)) * float(i)))
          y.append(np.cos((old_div(np.pi, 180.)) * float(i)))
      plt.plot(x, y, 'w-')
      x, y = [], []
      # the map will be a square of 1X1..this is how I erase the redundant area
      for j in range(1, 4):
          for i in range(0, 360):
              x.append(np.sin((old_div(np.pi, 180.)) * float(i))
                       * (1. + (old_div(float(j), 10.))))
              y.append(np.cos((old_div(np.pi, 180.)) * float(i))
                       * (1. + (old_div(float(j), 10.))))
          plt.plot(x, y, 'w-', linewidth=26)
          x, y = [], []
      # the axes
      plt.axis("equal")