def PintPars(datablock, araiblock, zijdblock, start, end, accept, **kwargs):
    """
     calculate the paleointensity magic parameters  make some definitions
    """
    if 'version' in list(kwargs.keys()) and kwargs['version'] == 3:
        meth_key = 'method_codes'
        beta_key = 'int_b_beta'
        temp_key, min_key, max_key = 'treat_temp', 'meas_step_min', 'meas_step_max'
        dc_theta_key, dc_phi_key = 'treat_dc_field_theta', 'treat_dc_field_phi'
        # convert dataframe to list of dictionaries
        datablock = datablock.to_dict('records')
        z_key = 'int_z'
        drats_key = 'int_drats'
        drat_key = 'int_drat'
        md_key = 'int_md'
        dec_key = 'dir_dec'
        inc_key = 'dir_inc'
        mad_key = 'int_mad_free'
        dang_key = 'int_dang'
        ptrm_key = 'int_n_ptrm'
        theta_key = 'int_theta'
        gamma_key = 'int_gamma'
        delta_key = 'int_delta'
        frac_key = 'int_frac'
        gmax_key = 'int_gmax'
        scat_key = 'int_scat'
    else:
        beta_key = 'specimen_b_beta'
        meth_key = 'magic_method_codes'
        temp_key, min_key, max_key = 'treatment_temp', 'measurement_step_min', 'measurement_step_max'
        z_key = 'specimen_z'
        drats_key = 'specimen_drats'
        drat_key = 'specimen_drat'
        md_key = 'specimen_md'
        dec_key = 'specimen_dec'
        inc_key = 'specimen_inc'
        mad_key = 'specimen_int_mad'
        dang_key = 'specimen_dang'
        ptrm_key = 'specimen_int_ptrm_n'
        theta_key = 'specimen_theta'
        gamma_key = 'specimen_gamma'
        delta_key = 'specimen_delta'
        frac_key = 'specimen_frac'
        gmax_key = 'specimen_gmax'
        scat_key = 'specimen_scat'

    first_Z, first_I, zptrm_check, ptrm_check, ptrm_tail = [], [], [], [], []
    methcode, ThetaChecks, DeltaChecks, GammaChecks = "", "", "", ""
    zptrm_check = []
    first_Z, first_I, ptrm_check, ptrm_tail, zptrm_check, GammaChecks = araiblock[
        0], araiblock[1], araiblock[2], araiblock[3], araiblock[4], araiblock[5]
    if len(araiblock) > 6:
        # used only for perpendicular method of paleointensity
        ThetaChecks = araiblock[6]
        # used only for perpendicular  method of paleointensity
        DeltaChecks = araiblock[7]
    xi, yi, diffcum = [], [], 0
    xiz, xzi, yiz, yzi = [], [], [], []
    Nptrm, dmax = 0, -1e-22
# check if even zero and infield steps
    if len(first_Z) > len(first_I):
        maxe = len(first_I) - 1
    else:
        maxe = len(first_Z) - 1
    if end == 0 or end > maxe:
        end = maxe
# get the MAD, DANG, etc. for directional data
    bstep = araiblock[0][start][0]
    estep = araiblock[0][end][0]
    zstart, zend = 0, len(zijdblock)
    for k in range(len(zijdblock)):
        zrec = zijdblock[k]
        if zrec[0] == bstep:
            zstart = k
        if zrec[0] == estep:
            zend = k
    PCA = domean(zijdblock, zstart, zend, 'DE-BFL')
    D, Diz, Dzi, Du = [], [], [], []  # list of NRM vectors, and separated by zi and iz
    for rec in zijdblock:
        D.append((rec[1], rec[2], rec[3]))
        Du.append((rec[1], rec[2]))
        if rec[4] == 1:
            Dzi.append((rec[1], rec[2]))  # if this is ZI step
        else:
            Diz.append((rec[1], rec[2]))  # if this is IZ step
# calculate the vector difference sum
    vds = dovds(D)
    b_zi, b_iz = [], []
# collect data included in ZigZag calculation
    if end + 1 >= len(first_Z):
        stop = end - 1
    else:
        stop = end
    for k in range(start, end + 1):
        for l in range(len(first_I)):
            irec = first_I[l]
            if irec[0] == first_Z[k][0]:
                xi.append(irec[3])
                yi.append(first_Z[k][3])
    pars, errcode = int_pars(xi, yi, vds)
    if errcode == 1:
        return pars, errcode
#    for k in range(start,end+1):
    for k in range(len(first_Z) - 1):
        for l in range(k):
            # only go down to 10% of NRM.....
            if old_div(first_Z[k][3], vds) > 0.1:
                irec = first_I[l]
                if irec[4] == 1 and first_I[l + 1][4] == 0:  # a ZI step
                    xzi = irec[3]
                    yzi = first_Z[k][3]
                    xiz = first_I[l + 1][3]
                    yiz = first_Z[k + 1][3]
                    slope = np.arctan2((yzi - yiz), (xiz - xzi))
                    r = np.sqrt((yzi - yiz)**2 + (xiz - xzi)**2)
                    if r > .1 * vds:
                        b_zi.append(slope)  # suppress noise
                elif irec[4] == 0 and first_I[l + 1][4] == 1:  # an IZ step
                    xiz = irec[3]
                    yiz = first_Z[k][3]
                    xzi = first_I[l + 1][3]
                    yzi = first_Z[k + 1][3]
                    slope = np.arctan2((yiz - yzi), (xzi - xiz))
                    r = np.sqrt((yiz - yzi)**2 + (xzi - xiz)**2)
                    if r > .1 * vds:
                        b_iz.append(slope)  # suppress noise
#
    ZigZag, Frat, Trat = -1, 0, 0
    if len(Diz) > 2 and len(Dzi) > 2:
        ZigZag = 0
        dizp = fisher_mean(Diz)  # get Fisher stats on IZ steps
        dzip = fisher_mean(Dzi)  # get Fisher stats on ZI steps
        dup = fisher_mean(Du)  # get Fisher stats on all steps
#
# if directions are TOO well grouped, can get false positive for ftest, so
# angles must be > 3 degrees apart.
#
        if angle([dizp['dec'], dizp['inc']], [dzip['dec'], dzip['inc']]) > 3.:
            F = (dup['n'] - 2.) * (dzip['r'] + dizp['r'] - dup['r']) / \
                (dup['n'] - dzip['r'] - dizp['r']
                 )  # Watson test for common mean
            nf = 2. * (dup['n'] - 2.)  # number of degees of freedom
            ftest = fcalc(2, nf)
            Frat = old_div(F, ftest)
            if Frat > 1.:
                ZigZag = Frat  # fails zigzag on directions
                methcode = "SM-FTEST"
# now do slopes
    if len(b_zi) > 2 and len(b_iz) > 2:
        bzi_m, bzi_sig = gausspars(b_zi)  # mean, std dev
        biz_m, biz_sig = gausspars(b_iz)
        n_zi = float(len(b_zi))
        n_iz = float(len(b_iz))
        b_diff = abs(bzi_m - biz_m)  # difference in means
#
# avoid false positives - set 3 degree slope difference here too
        if b_diff > 3 * np.pi / 180.:
            nf = n_zi + n_iz - 2.  # degrees of freedom
            svar = old_div(((n_zi - 1.) * bzi_sig**2 +
                            (n_iz - 1.) * biz_sig**2), nf)
            T = old_div((b_diff), np.sqrt(
                svar * (old_div(1.0, n_zi) + old_div(1.0, n_iz))))  # student's t
            ttest = tcalc(nf, .05)  # t-test at 95% conf.
            Trat = old_div(T, ttest)
            if Trat > 1 and Trat > Frat:
                ZigZag = Trat  # fails zigzag on directions
                methcode = "SM-TTEST"
    pars[z_key] = ZigZag
    pars[meth_key] = methcode
# do drats
    if len(ptrm_check) != 0:
        diffcum, drat_max = 0, 0
        for prec in ptrm_check:
            step = prec[0]
            endbak = end
            zend = end
            while zend > len(zijdblock) - 1:
                zend = zend - 2  # don't count alteration that happens after this step
            if step < zijdblock[zend][0]:
                Nptrm += 1
                for irec in first_I:
                    if irec[0] == step:
                        break
                diffcum += prec[3] - irec[3]
                if abs(prec[3] - irec[3]) > drat_max:
                    drat_max = abs(prec[3] - irec[3])
        pars[drats_key] = (100 * abs(diffcum) / first_I[zend][3])
        pars[drat_key] = (100 * abs(drat_max) / first_I[zend][3])
    elif len(zptrm_check) != 0:
        diffcum = 0
        for prec in zptrm_check:
            step = prec[0]
            endbak = end
            zend = end
            while zend > len(zijdblock) - 1:
                zend = zend - 1
            if step < zijdblock[zend][0]:
                Nptrm += 1
                for irec in first_I:
                    if irec[0] == step:
                        break
                diffcum += prec[3] - irec[3]
        pars[drats_key] = (100 * abs(diffcum) / first_I[zend][3])
    else:
        pars[drats_key] = -1
        pars[drat_key] = -1
# and the pTRM tails
    if len(ptrm_tail) != 0:
        for trec in ptrm_tail:
            step = trec[0]
            for irec in first_I:
                if irec[0] == step:
                    break
            if abs(trec[3]) > dmax:
                dmax = abs(trec[3])
        pars[md_key] = (100 * dmax / vds)
    else:
        pars[md_key] = -1
    pars[min_key] = bstep
    pars[max_key] = estep
    pars[dec_key] = PCA["specimen_dec"]
    pars[inc_key] = PCA["specimen_inc"]
    pars[mad_key] = PCA["specimen_mad"]
    pars[dang_key] = PCA["specimen_dang"]
    pars[ptrm_key] = Nptrm
# and the ThetaChecks
    if ThetaChecks != "":
        t = 0
        for theta in ThetaChecks:
            if theta[0] >= bstep and theta[0] <= estep and theta[1] > t:
                t = theta[1]
        pars[theta_key] = t
    else:
        pars[theta_key] = -1
# and the DeltaChecks
    if DeltaChecks != "":
        d = 0
        for delta in DeltaChecks:
            if delta[0] >= bstep and delta[0] <= estep and delta[1] > d:
                d = delta[1]
        pars[delta_key]
    else:
        pars[delta_key] = -1
    pars[gamma_key] = -1
    if GammaChecks != "":
        for gamma in GammaChecks:
            if gamma[0] <= estep:
                pars['specimen_gamma'] = gamma[1]

    # --------------------------------------------------------------
    # From here added By Ron Shaar 11-Dec 2012
    # New parameters defined in Shaar and Tauxe (2012):
    # FRAC (specimen_frac) - ranges from 0. to 1.
    # SCAT (specimen_scat) - takes 1/0
    # gap_max (specimen_gmax) - ranges from 0. to 1.
    # --------------------------------------------------------------

    # --------------------------------------------------------------
    # FRAC is similar to Fvds, but the numerator is the vds fraction:
    # FRAC= [ vds (start,end)] / total vds ]
    # gap_max= max [ (vector difference) /  vds (start,end)]
    # --------------------------------------------------------------

    # collect all zijderveld data to arrays and calculate VDS

    z_temperatures = [row[0] for row in zijdblock]
    zdata = []                # array of zero-fields measurements in Cartezian coordinates
    # array of vector differences (for vds calculation)
    vector_diffs = []
    NRM = zijdblock[0][3]     # NRM

    for k in range(len(zijdblock)):
        DIR = [zijdblock[k][1], zijdblock[k][2], old_div(zijdblock[k][3], NRM)]
        cart = dir2cart(DIR)
        zdata.append(np.array([cart[0], cart[1], cart[2]]))
        if k > 0:
            vector_diffs.append(
                np.sqrt(sum((np.array(zdata[-2]) - np.array(zdata[-1]))**2)))
    # last vector difference: from the last point to the origin.
    vector_diffs.append(np.sqrt(sum(np.array(zdata[-1])**2)))
    vds = sum(vector_diffs)  # vds calculation
    zdata = np.array(zdata)
    vector_diffs = np.array(vector_diffs)

    # calculate the vds within the chosen segment
    vector_diffs_segment = vector_diffs[zstart:zend]
    # FRAC calculation
    FRAC = old_div(sum(vector_diffs_segment), vds)
    pars[frac_key] = FRAC

    # gap_max calculation
    max_FRAC_gap = max(
        old_div(vector_diffs_segment, sum(vector_diffs_segment)))
    pars[gmax_key] = max_FRAC_gap

    # ---------------------------------------------------------------------
    # Calculate the "scat box"
    # all data-points, pTRM checks, and tail-checks, should be inside a "scat box"
    # ---------------------------------------------------------------------

    # intialization
    # fail scat due to arai plot data points
    pars["fail_arai_beta_box_scatter"] = False
    pars["fail_ptrm_beta_box_scatter"] = False  # fail scat due to pTRM checks
    pars["fail_tail_beta_box_scatter"] = False  # fail scat due to tail checks
    pars[scat_key] = "t"  # Pass by default

    # --------------------------------------------------------------
    # collect all Arai plot data points in arrays

    x_Arai, y_Arai, t_Arai, steps_Arai = [], [], [], []
    NRMs = araiblock[0]
    PTRMs = araiblock[1]
    ptrm_checks = araiblock[2]
    ptrm_tail = araiblock[3]

    PTRMs_temperatures = [row[0] for row in PTRMs]
    NRMs_temperatures = [row[0] for row in NRMs]
    NRM = NRMs[0][3]

    for k in range(len(NRMs)):
        index_pTRMs = PTRMs_temperatures.index(NRMs[k][0])
        x_Arai.append(old_div(PTRMs[index_pTRMs][3], NRM))
        y_Arai.append(old_div(NRMs[k][3], NRM))
        t_Arai.append(NRMs[k][0])
        if NRMs[k][4] == 1:
            steps_Arai.append('ZI')
        else:
            steps_Arai.append('IZ')
    x_Arai = np.array(x_Arai)
    y_Arai = np.array(y_Arai)

    # --------------------------------------------------------------
    # collect all pTRM check to arrays

    x_ptrm_check, y_ptrm_check, ptrm_checks_temperatures, = [], [], []
    x_ptrm_check_starting_point, y_ptrm_check_starting_point, ptrm_checks_starting_temperatures = [], [], []

    for k in range(len(ptrm_checks)):
        if ptrm_checks[k][0] in NRMs_temperatures:
            # find the starting point of the pTRM check:
            for i in range(len(datablock)):
                rec = datablock[i]
                if "LT-PTRM-I" in rec[meth_key] and float(rec[temp_key]) == ptrm_checks[k][0]:
                    starting_temperature = (float(datablock[i - 1][temp_key]))
                    try:
                        index = t_Arai.index(starting_temperature)
                        x_ptrm_check_starting_point.append(x_Arai[index])
                        y_ptrm_check_starting_point.append(y_Arai[index])
                        ptrm_checks_starting_temperatures.append(
                            starting_temperature)

                        index_zerofield = zerofield_temperatures.index(
                            ptrm_checks[k][0])
                        x_ptrm_check.append(old_div(ptrm_checks[k][3], NRM))
                        y_ptrm_check.append(
                            old_div(zerofields[index_zerofield][3], NRM))
                        ptrm_checks_temperatures.append(ptrm_checks[k][0])

                        break
                    except:
                        pass

    x_ptrm_check_starting_point = np.array(x_ptrm_check_starting_point)
    y_ptrm_check_starting_point = np.array(y_ptrm_check_starting_point)
    ptrm_checks_starting_temperatures = np.array(
        ptrm_checks_starting_temperatures)
    x_ptrm_check = np.array(x_ptrm_check)
    y_ptrm_check = np.array(y_ptrm_check)
    ptrm_checks_temperatures = np.array(ptrm_checks_temperatures)

    # --------------------------------------------------------------
    # collect tail checks to arrays

    x_tail_check, y_tail_check, tail_check_temperatures = [], [], []
    x_tail_check_starting_point, y_tail_check_starting_point, tail_checks_starting_temperatures = [], [], []

    for k in range(len(ptrm_tail)):
        if ptrm_tail[k][0] in NRMs_temperatures:

            # find the starting point of the pTRM check:
            for i in range(len(datablock)):
                rec = datablock[i]
                if "LT-PTRM-MD" in rec[meth_key] and float(rec[temp_key]) == ptrm_tail[k][0]:
                    starting_temperature = (float(datablock[i - 1][temp_key]))
                    try:

                        index = t_Arai.index(starting_temperature)
                        x_tail_check_starting_point.append(x_Arai[index])
                        y_tail_check_starting_point.append(y_Arai[index])
                        tail_checks_starting_temperatures.append(
                            starting_temperature)

                        index_infield = infield_temperatures.index(
                            ptrm_tail[k][0])
                        x_tail_check.append(
                            old_div(infields[index_infield][3], NRM))
                        y_tail_check.append(
                            old_div(ptrm_tail[k][3], NRM) + old_div(zerofields[index_infield][3], NRM))
                        tail_check_temperatures.append(ptrm_tail[k][0])

                        break
                    except:
                        pass

    x_tail_check = np.array(x_tail_check)
    y_tail_check = np.array(y_tail_check)
    tail_check_temperatures = np.array(tail_check_temperatures)
    x_tail_check_starting_point = np.array(x_tail_check_starting_point)
    y_tail_check_starting_point = np.array(y_tail_check_starting_point)
    tail_checks_starting_temperatures = np.array(
        tail_checks_starting_temperatures)

    # --------------------------------------------------------------
    # collect the chosen segment in the Arai plot to arrays

    x_Arai_segment = x_Arai[start:end + 1]  # chosen segent in the Arai plot
    y_Arai_segment = y_Arai[start:end + 1]  # chosen segent in the Arai plot

    # --------------------------------------------------------------
    # collect pTRM checks in segment to arrays
    # notice, this is different than the conventional DRATS.
    # for scat calculation we take only the pTRM checks which were carried out
    # before reaching the highest temperature in the chosen segment

    x_ptrm_check_for_SCAT, y_ptrm_check_for_SCAT = [], []
    for k in range(len(ptrm_checks_temperatures)):
        if ptrm_checks_temperatures[k] >= pars[min_key] and ptrm_checks_starting_temperatures <= pars[max_key]:
            x_ptrm_check_for_SCAT.append(x_ptrm_check[k])
            y_ptrm_check_for_SCAT.append(y_ptrm_check[k])

    x_ptrm_check_for_SCAT = np.array(x_ptrm_check_for_SCAT)
    y_ptrm_check_for_SCAT = np.array(y_ptrm_check_for_SCAT)

    # --------------------------------------------------------------
    # collect Tail checks in segment to arrays
    # for scat calculation we take only the tail checks which were carried out
    # before reaching the highest temperature in the chosen segment

    x_tail_check_for_SCAT, y_tail_check_for_SCAT = [], []

    for k in range(len(tail_check_temperatures)):
        if tail_check_temperatures[k] >= pars[min_key] and tail_checks_starting_temperatures[k] <= pars[max_key]:
            x_tail_check_for_SCAT.append(x_tail_check[k])
            y_tail_check_for_SCAT.append(y_tail_check[k])

    x_tail_check_for_SCAT = np.array(x_tail_check_for_SCAT)
    y_tail_check_for_SCAT = np.array(y_tail_check_for_SCAT)

    # --------------------------------------------------------------
    # calculate the lines that define the scat box:

    # if threshold value for beta is not defined, then scat cannot be calculated (pass)
    # in this case, scat pass
    if beta_key in list(accept.keys()) and accept[beta_key] != "":
        b_beta_threshold = float(accept[beta_key])
        b = pars[b_key]             # best fit line
        cm_x = np.mean(np.array(x_Arai_segment))  # x center of mass
        cm_y = np.mean(np.array(y_Arai_segment))  # y center of mass
        a = cm_y - b * cm_x

        # lines with slope = slope +/- 2*(specimen_b_beta)

        two_sigma_beta_threshold = 2 * b_beta_threshold
        two_sigma_slope_threshold = abs(two_sigma_beta_threshold * b)

        # a line with a  shallower  slope  (b + 2*beta*b) passing through the center of mass
        # y=a1+b1x
        b1 = b + two_sigma_slope_threshold
        a1 = cm_y - b1 * cm_x

        # bounding line with steeper  slope (b - 2*beta*b) passing through the center of mass
        # y=a2+b2x
        b2 = b - two_sigma_slope_threshold
        a2 = cm_y - b2 * cm_x

        # lower bounding line of the 'beta box'
        # y=intercept1+slop1x
        slop1 = old_div(a1, ((old_div(a2, b2))))
        intercept1 = a1

        # higher bounding line of the 'beta box'
        # y=intercept2+slop2x

        slop2 = old_div(a2, ((old_div(a1, b1))))
        intercept2 = a2

        pars['specimen_scat_bounding_line_high'] = [intercept2, slop2]
        pars['specimen_scat_bounding_line_low'] = [intercept1, slop1]

        # --------------------------------------------------------------
        # check if the Arai data points are in the 'box'

        # the two bounding lines
        ymin = intercept1 + x_Arai_segment * slop1
        ymax = intercept2 + x_Arai_segment * slop2

        # arrays of "True" or "False"
        check_1 = y_Arai_segment > ymax
        check_2 = y_Arai_segment < ymin

        # check if at least one "True"
        if (sum(check_1) + sum(check_2)) > 0:
            pars["fail_arai_beta_box_scatter"] = True

        # --------------------------------------------------------------
        # check if the pTRM checks data points are in the 'box'

        if len(x_ptrm_check_for_SCAT) > 0:

            # the two bounding lines
            ymin = intercept1 + x_ptrm_check_for_SCAT * slop1
            ymax = intercept2 + x_ptrm_check_for_SCAT * slop2

            # arrays of "True" or "False"
            check_1 = y_ptrm_check_for_SCAT > ymax
            check_2 = y_ptrm_check_for_SCAT < ymin

            # check if at least one "True"
            if (sum(check_1) + sum(check_2)) > 0:
                pars["fail_ptrm_beta_box_scatter"] = True

        # --------------------------------------------------------------
        # check if the tail checks data points are in the 'box'

        if len(x_tail_check_for_SCAT) > 0:

            # the two bounding lines
            ymin = intercept1 + x_tail_check_for_SCAT * slop1
            ymax = intercept2 + x_tail_check_for_SCAT * slop2

            # arrays of "True" or "False"
            check_1 = y_tail_check_for_SCAT > ymax
            check_2 = y_tail_check_for_SCAT < ymin

            # check if at least one "True"
            if (sum(check_1) + sum(check_2)) > 0:
                pars["fail_tail_beta_box_scatter"] = True

        # --------------------------------------------------------------
        # check if specimen_scat is PASS or FAIL:

        if pars["fail_tail_beta_box_scatter"] or pars["fail_ptrm_beta_box_scatter"] or pars["fail_arai_beta_box_scatter"]:
            pars[scat_key] = 'f'
        else:
            pars[scat_key] = 't'

    return pars, 0