docs: def plot_dir(ZED, pars, datablock, angle):
      """
      function to put the great circle on the equal area projection
      and plot start and end points of calculation
  
      DEPRECATED (used in zeq_magic)
      """
  #
  # find start and end points from datablock
  #
      if pars["calculation_type"] == 'DE-FM':
          x, y = [], []
          plt.figure(num=ZED['eqarea'])
          XY = pmag.dimap(pars["specimen_dec"], pars["specimen_inc"])
          x.append(XY[0])
          y.append(XY[1])
          plt.scatter(x, y, marker='^', s=80, c='r')
          plt.show()
          return
      StartDir, EndDir = [0, 0, 1.], [0, 0, 1.]
      for rec in datablock:
          if rec[0] == pars["measurement_step_min"]:
              StartDir[0] = rec[1]
              StartDir[1] = rec[2]
              if pars["specimen_direction_type"] == 'l':
                  StartDir[2] = rec[3]/datablock[0][3]
          if rec[0] == pars["measurement_step_max"]:
              EndDir[0] = rec[1]
              EndDir[1] = rec[2]
              if pars["specimen_direction_type"] == 'l':
                  EndDir[2] = rec[3]/datablock[0][3]
  
  #
  #  put them on the plots
  #
      x, y, z, pole = [], [], [], []
      if pars["calculation_type"] != 'DE-BFP':
          plt.figure(num=ZED['eqarea'])
          XY = pmag.dimap(pars["specimen_dec"], pars["specimen_inc"])
          x.append(XY[0])
          y.append(XY[1])
          plt.scatter(x, y, marker='d', s=80, c='b')
          x, y, z = [], [], []
          StartDir[0] = StartDir[0] - angle
          EndDir[0] = EndDir[0] - angle
          XYZs = pmag.dir2cart(StartDir)
          x.append(XYZs[0])
          y.append(XYZs[1])
          z.append(XYZs[2])
          XYZe = pmag.dir2cart(EndDir)
          x.append(XYZe[0])
          y.append(XYZe[1])
          z.append(XYZe[2])
          plt.figure(num=ZED['zijd'])
          plt.scatter(x, y, marker='d', s=80, c='g')
          plt.scatter(x, z, marker='d', s=80, c='g')
          plt.scatter(x, y, marker='o', c='r', s=20)
          plt.scatter(x, z, marker='s', c='w', s=20)
  #
  # put on best fit line
  # new way (from Jeff Gee's favorite website http://GET THIS):
  #      P1=pmag.dir2cart([(pars["specimen_dec"]-angle),pars["specimen_inc"],1.]) #  princ comp.
  #      P2=pmag.dir2cart([(pars["specimen_dec"]-angle-180.),-pars["specimen_inc"],1.]) # antipode of princ comp.
  #      P21,Ps,Pe,Xs,Xe=[],[],[],[],[]
  #      for i in range(3):
  #          P21.append(P2[i]-P1[i])
  #          Ps.append(XYZs[i]-P1[i])
  #          Pe.append(XYZe[i]-P1[i])
  #      norm=pmag.cart2dir(P21)[2]
  #      us=(Ps[0]*P21[0]+Ps[1]*P21[1]+Ps[2]*P21[2])/(norm**2)
  #      ue=(Pe[0]*P21[0]+Pe[1]*P21[1]+Pe[2]*P21[2])/(norm**2)
  #      px,py,pz=[],[],[]
  #      for i in range(3):
  #          Xs.append(P1[i]+us*(P2[i]-P1[i]))
  #          Xe.append(P1[i]+ue*(P2[i]-P1[i]))
  #   old way:
          cm = pars["center_of_mass"]
          if cm != [0., 0., 0.]:
              cm = np.array(pars["center_of_mass"])/datablock[0][3]
              cmDir = pmag.cart2dir(cm)
              cmDir[0] = cmDir[0] - angle
              cm = pmag.dir2cart(cmDir)
              diff = []
              for i in range(3):
                  diff.append(XYZe[i] - XYZs[i])
              R = np.sqrt(diff[0]**2 + diff[1]**2 + diff[2]**2)
              P = pmag.dir2cart(
                  ((pars["specimen_dec"] - angle), pars["specimen_inc"], R/2.5))
              px, py, pz = [], [], []
              px.append((cm[0] + P[0]))
              py.append((cm[1] + P[1]))
              pz.append((cm[2] + P[2]))
              px.append((cm[0] - P[0]))
              py.append((cm[1] - P[1]))
              pz.append((cm[2] - P[2]))
  
          plt.plot(px, py, 'g', linewidth=2)
          plt.plot(px, pz, 'g', linewidth=2)
          plt.axis("equal")
      else:
          plt.figure(num=ZED['eqarea'])
          XY = pmag.dimap(StartDir[0], StartDir[1])
          x.append(XY[0])
          y.append(XY[1])
          XY = pmag.dimap(EndDir[0], EndDir[1])
          x.append(XY[0])
          y.append(XY[1])
          plt.scatter(x, y, marker='d', s=80, c='b')
          pole.append(pars["specimen_dec"])
          pole.append(pars["specimen_inc"])
          plot_circ(ZED['eqarea'], pole, 90., 'g')
          plt.xlim((-1., 1.))
          plt.ylim((-1., 1.))
          plt.axis("equal")