def domean(data, start, end, calculation_type):
    """
    Gets average direction using Fisher or principal component analysis (line
    or plane) methods

    Parameters
    ----------
    data : nest list of data: [[treatment,dec,inc,int,quality],...]
    start : step being used as start of fit (often temperature minimum)
    end : step being used as end of fit (often temperature maximum)
    calculation_type : string describing type of calculation to be made
    'DE-BFL' (line), 'DE-BFL-A' (line-anchored), 'DE-BFL-O' (line-with-origin),
    'DE-BFP' (plane), 'DE-FM' (Fisher mean)

    Returns
    -------
    mpars : dictionary with the keys "specimen_n","measurement_step_min",
    "measurement_step_max","specimen_mad","specimen_dec","specimen_inc"
    """
    mpars = {}
    datablock = []
    start0, end0 = start, end
    # indata = [rec.append('g') if len(rec)<6 else rec for rec in indata] #
    # this statement doesn't work!
    indata = []
    for rec in data:
        if len(rec) < 6:
            rec.append('g')
        indata.append(rec)
    if indata[start0][5] == 'b':
        print("Can't select 'bad' point as start for PCA")
    flags = [x[5] for x in indata]
    bad_before_start = flags[:start0].count('b')
    bad_in_mean = flags[start0:end0 + 1].count('b')
    start = start0 - bad_before_start
    end = end0 - bad_before_start - bad_in_mean
    datablock = [x for x in indata if x[5] == 'g']
    if indata[start0] != datablock[start]:
        print('problem removing bad data in pmag.domean start of datablock shifted:\norigional: %d\nafter removal: %d' % (
            start0, indata.index(datablock[start])))
    if indata[end0] != datablock[end]:
        print('problem removing bad data in pmag.domean end of datablock shifted:\norigional: %d\nafter removal: %d' % (
            end0, indata.index(datablock[end])))
    mpars["calculation_type"] = calculation_type
    rad = old_div(np.pi, 180.)
    if end > len(datablock) - 1 or end < start:
        end = len(datablock) - 1
    control, data, X, Nrec = [], [], [], float(end - start + 1)
    cm = [0., 0., 0.]
#
#  get cartesian coordinates
#
    fdata = []
    for k in range(start, end + 1):
        if calculation_type == 'DE-BFL' or calculation_type == 'DE-BFL-A' or calculation_type == 'DE-BFL-O':  # best-fit line
            data = [datablock[k][1], datablock[k][2], datablock[k][3]]
        else:
            data = [datablock[k][1], datablock[k][2], 1.0]  # unit weight
        fdata.append(data)
        cart = dir2cart(data)
        X.append(cart)
    if calculation_type == 'DE-BFL-O':  # include origin as point
        X.append([0., 0., 0.])
        # pass
    if calculation_type == 'DE-FM':  # for fisher means
        fpars = fisher_mean(fdata)
        mpars["specimen_direction_type"] = 'l'
        mpars["specimen_dec"] = fpars["dec"]
        mpars["specimen_inc"] = fpars["inc"]
        mpars["specimen_alpha95"] = fpars["alpha95"]
        mpars["specimen_n"] = fpars["n"]
        mpars["specimen_r"] = fpars["r"]
        mpars["measurement_step_min"] = indata[start0][0]
        mpars["measurement_step_max"] = indata[end0][0]
        mpars["center_of_mass"] = cm
        mpars["specimen_dang"] = -1
        return mpars
#
#   get center of mass for principal components (DE-BFL or DE-BFP)
#
    for cart in X:
        for l in range(3):
            cm[l] += old_div(cart[l], Nrec)
    mpars["center_of_mass"] = cm

#
#   transform to center of mass (if best-fit line)
#
    if calculation_type != 'DE-BFP':
        mpars["specimen_direction_type"] = 'l'
    if calculation_type == 'DE-BFL' or calculation_type == 'DE-BFL-O':  # not for planes or anchored lines
        for k in range(len(X)):
            for l in range(3):
                X[k][l] = X[k][l] - cm[l]
    else:
        mpars["specimen_direction_type"] = 'p'

#
#   put in T matrix
#
    T = np.array(Tmatrix(X))
#
#   get sorted evals/evects
#
    t, V = tauV(T)
    if t == []:
        mpars["specimen_direction_type"] = "Error"
        print("Error in calculation")
        return mpars
    v1, v3 = V[0], V[2]
    if t[2] < 0:
        t[2] = 0  # make positive
    if calculation_type == 'DE-BFL-A':
        Dir, R = vector_mean(fdata)
        mpars["specimen_direction_type"] = 'l'
        mpars["specimen_dec"] = Dir[0]
        mpars["specimen_inc"] = Dir[1]
        mpars["specimen_n"] = len(fdata)
        mpars["measurement_step_min"] = indata[start0][0]
        mpars["measurement_step_max"] = indata[end0][0]
        mpars["center_of_mass"] = cm
        s1 = np.sqrt(t[0])
        MAD = old_div(np.arctan(old_div(np.sqrt(t[1] + t[2]), s1)), rad)
        if np.iscomplexobj(MAD):
            MAD = MAD.real
        # I think this is how it is done - i never anchor the "PCA" - check
        mpars["specimen_mad"] = MAD
        return mpars
    if calculation_type != 'DE-BFP':
        #
        #   get control vector for principal component direction
        #
        rec = [datablock[start][1], datablock[start][2], datablock[start][3]]
        P1 = dir2cart(rec)
        rec = [datablock[end][1], datablock[end][2], datablock[end][3]]
        P2 = dir2cart(rec)
#
#   get right direction along principal component
##
        for k in range(3):
            control.append(P1[k] - P2[k])
        # changed by rshaar
        # control is taken as the center of mass
        # control=cm

        dot = 0
        for k in range(3):
            dot += v1[k] * control[k]
        if dot < -1:
            dot = -1
        if dot > 1:
            dot = 1
        if np.arccos(dot) > old_div(np.pi, 2.):
            for k in range(3):
                v1[k] = -v1[k]
#   get right direction along principal component
#
        s1 = np.sqrt(t[0])
        Dir = cart2dir(v1)
        MAD = old_div(np.arctan(old_div(np.sqrt(t[1] + t[2]), s1)), rad)
        if np.iscomplexobj(MAD):
            MAD = MAD.real
    if calculation_type == "DE-BFP":
        Dir = cart2dir(v3)
        MAD = old_div(
            np.arctan(np.sqrt(old_div(t[2], t[1]) + old_div(t[2], t[0]))), rad)
        if np.iscomplexobj(MAD):
            MAD = MAD.real
#
#   get angle with  center of mass
#
    CMdir = cart2dir(cm)
    Dirp = [Dir[0], Dir[1], 1.]
    dang = angle(CMdir, Dirp)
    mpars["specimen_dec"] = Dir[0]
    mpars["specimen_inc"] = Dir[1]
    mpars["specimen_mad"] = MAD
    # mpars["specimen_n"]=int(Nrec)
    mpars["specimen_n"] = len(X)
    mpars["specimen_dang"] = dang[0]
    mpars["measurement_step_min"] = indata[start0][0]
    mpars["measurement_step_max"] = indata[end0][0]
    return mpars
Gets average direction using Fisher or principal component analysis (line
    or plane) methods

    Parameters
    ----------
    data : nest list of data: [[treatment,dec,inc,int,quality],...]
    start : step being used as start of fit (often temperature minimum)
    end : step being used as end of fit (often temperature maximum)
    calculation_type : string describing type of calculation to be made
    'DE-BFL' (line), 'DE-BFL-A' (line-anchored), 'DE-BFL-O' (line-with-origin),
    'DE-BFP' (plane), 'DE-FM' (Fisher mean)

    Returns
    -------
    mpars : dictionary with the keys "specimen_n","measurement_step_min",
    "measurement_step_max","specimen_mad","specimen_dec","specimen_inc"