docs: def find_f(data):
      """
      Given a distribution of directions, this function determines parameters
      (elongation, inclination, flattening factor, and elongation direction) that
      are consistent with the TK03 secular variation model.
  
      Parameters
      ----------
      data : array of declination, inclination pairs
          (e.g. np.array([[140,21],[127,23],[142,19],[136,22]]))
  
      Returns
      ---------
      Es : list of elongation values
      Is : list of inclination values
      Fs : list of flattening factors
      V2s : list of elongation directions (relative to the distribution)
  
      The function will return a zero list ([0]) for each of these parameters if the directions constitute a pathological distribution.
  
      Examples
      ---------
      >>> directions = np.array([[140,21],[127,23],[142,19],[136,22]])
      >>> Es, Is, Fs, V2s = pmag.find_f(directions)
      """
      rad = np.pi/180.
      Es, Is, Fs, V2s = [], [], [], []
      ppars = doprinc(data)
      D = ppars['dec']
      Decs, Incs = data.transpose()[0], data.transpose()[1]
      Tan_Incs = np.tan(Incs * rad)
      for f in np.arange(1., .2, -.01):
          U = old_div(np.arctan((old_div(1., f)) * Tan_Incs), rad)
          fdata = np.array([Decs, U]).transpose()
          ppars = doprinc(fdata)
          Fs.append(f)
          Es.append(old_div(ppars["tau2"], ppars["tau3"]))
          ang = angle([D, 0], [ppars["V2dec"], 0])
          if 180. - ang < ang:
              ang = 180. - ang
          V2s.append(ang)
          Is.append(abs(ppars["inc"]))
          if EI(abs(ppars["inc"])) <= Es[-1]:
              del Es[-1]
              del Is[-1]
              del Fs[-1]
              del V2s[-1]
              if len(Fs) > 0:
                  for f in np.arange(Fs[-1], .2, -.005):
                      U = old_div(np.arctan((old_div(1., f)) * Tan_Incs), rad)
                      fdata = np.array([Decs, U]).transpose()
                      ppars = doprinc(fdata)
                      Fs.append(f)
                      Es.append(old_div(ppars["tau2"], ppars["tau3"]))
                      Is.append(abs(ppars["inc"]))
                      ang = angle([D, 0], [ppars["V2dec"], 0])
                      if 180. - ang < ang:
                          ang = 180. - ang
                      V2s.append(ang)
                      if EI(abs(ppars["inc"])) <= Es[-1]:
                          return Es, Is, Fs, V2s
      return [0], [0], [0], [0]