docs: def doincfish(inc):
      """
      gets fisher mean inc from inc only data
      input: list of inclination values
      output: dictionary of
          'n' : number of inclination values supplied
          'ginc' : gaussian mean of inclinations
          'inc' : estimated Fisher mean
          'r' : estimated Fisher R value
          'k' : estimated Fisher kappa
          'alpha95' : estimated fisher alpha_95
          'csd' : estimated circular standard deviation
      """
      rad, SCOi, SSOi = old_div(np.pi, 180.), 0., 0.  # some definitions
      abinc = []
      for i in inc:
          abinc.append(abs(i))
      MI, std = gausspars(abinc)  # get mean inc and standard deviation
      fpars = {}
      N = len(inc)  # number of data
      fpars['n'] = N
      fpars['ginc'] = MI
      if MI < 30:
          fpars['inc'] = MI
          fpars['k'] = 0
          fpars['alpha95'] = 0
          fpars['csd'] = 0
          fpars['r'] = 0
          print('WARNING: mean inc < 30, returning gaussian mean')
          return fpars
      for i in inc:  # sum over all incs (but take only positive inc)
          coinc = (90. - abs(i)) * rad
          SCOi += np.cos(coinc)
          SSOi += np.sin(coinc)
      Oo = (90.0 - MI) * rad  # first guess at mean
      SCFlag = -1  # sign change flag
      epsilon = float(N) * np.cos(Oo)  # RHS of zero equations
      epsilon += (np.sin(Oo)**2 - np.cos(Oo)**2) * SCOi
      epsilon -= 2. * np.sin(Oo) * np.cos(Oo) * SSOi
      while SCFlag < 0:  # loop until cross zero
          if MI > 0:
              Oo -= (.01 * rad)  # get steeper
          if MI < 0:
              Oo += (.01 * rad)  # get shallower
          prev = epsilon
          epsilon = float(N) * np.cos(Oo)  # RHS of zero equations
          epsilon += (np.sin(Oo)**2. - np.cos(Oo)**2.) * SCOi
          epsilon -= 2. * np.sin(Oo) * np.cos(Oo) * SSOi
          if abs(epsilon) > abs(prev):
              MI = -1 * MI  # reverse direction
          if epsilon * prev < 0:
              SCFlag = 1  # changed sign
      S, C = 0., 0.  # initialize for summation
      for i in inc:
          coinc = (90. - abs(i)) * rad
          S += np.sin(Oo - coinc)
          C += np.cos(Oo - coinc)
      k = old_div((N - 1.), (2. * (N - C)))
      Imle = 90. - (old_div(Oo, rad))
      fpars["inc"] = Imle
      fpars["r"], R = 2. * C - N, 2 * C - N
      fpars["k"] = k
      f = fcalc(2, N - 1)
      a95 = 1. - (0.5) * (old_div(S, C))**2 - (old_div(f, (2. * C * k)))
  #    b=20.**(1./(N-1.)) -1.
  #    a=1.-b*(N-R)/R
  #    a95=np.arccos(a)*180./np.pi
      csd = old_div(81., np.sqrt(k))
      fpars["alpha95"] = a95
      fpars["csd"] = csd
      return fpars