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
def PintPars ( datablock , araiblock , zijdblock , start , end , accept , ** kwargs ) : 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' 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 : ThetaChecks = araiblock [ 6 ] DeltaChecks = araiblock [ 7 ] xi , yi , diffcum = [ ] , [ ] , 0 xiz , xzi , yiz , yzi = [ ] , [ ] , [ ] , [ ] Nptrm , dmax = 0 , - 1e-22 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 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 = [ ] , [ ] , [ ] , [ ] 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 ] ) ) else : Diz . append ( ( rec [ 1 ] , rec [ 2 ] ) ) vds = dovds ( D ) b_zi , b_iz = [ ] , [ ] 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 ( len ( first_Z ) - 1 ) : for l in range ( k ) : 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 : 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 ) elif irec [ 4 ] == 0 and first_I [ l + 1 ] [ 4 ] == 1 : 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 ) ZigZag , Frat , Trat = - 1 , 0 , 0 if len ( Diz ) > 2 and len ( Dzi ) > 2 : ZigZag = 0 dizp = fisher_mean ( Diz ) dzip = fisher_mean ( Dzi ) dup = fisher_mean ( Du ) 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' ] ) nf = 2. * ( dup [ 'n' ] - 2. ) ftest = fcalc ( 2 , nf ) Frat = old_div ( F , ftest ) if Frat > 1. : ZigZag = Frat methcode = "SM-FTEST" if len ( b_zi ) > 2 and len ( b_iz ) > 2 : bzi_m , bzi_sig = gausspars ( b_zi ) 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 ) if b_diff > 3 * np . pi / 180. : nf = n_zi + n_iz - 2. 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 ) ) ) ) ttest = tcalc ( nf , .05 ) Trat = old_div ( T , ttest ) if Trat > 1 and Trat > Frat : ZigZag = Trat methcode = "SM-TTEST" pars [ z_key ] = ZigZag pars [ meth_key ] = methcode 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 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 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 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 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 ] z_temperatures = [ row [ 0 ] for row in zijdblock ] zdata = [ ] vector_diffs = [ ] NRM = zijdblock [ 0 ] [ 3 ] 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 ) ) ) vector_diffs . append ( np . sqrt ( sum ( np . array ( zdata [ - 1 ] ) ** 2 ) ) ) vds = sum ( vector_diffs ) zdata = np . array ( zdata ) vector_diffs = np . array ( vector_diffs ) vector_diffs_segment = vector_diffs [ zstart : zend ] FRAC = old_div ( sum ( vector_diffs_segment ) , vds ) pars [ frac_key ] = FRAC max_FRAC_gap = max ( old_div ( vector_diffs_segment , sum ( vector_diffs_segment ) ) ) pars [ gmax_key ] = max_FRAC_gap pars [ "fail_arai_beta_box_scatter" ] = False pars [ "fail_ptrm_beta_box_scatter" ] = False pars [ "fail_tail_beta_box_scatter" ] = False pars [ scat_key ] = "t" 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 ) 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 : 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 ) 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 : 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 ) x_Arai_segment = x_Arai [ start : end + 1 ] y_Arai_segment = y_Arai [ start : end + 1 ] 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 ) 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 ) if beta_key in list ( accept . keys ( ) ) and accept [ beta_key ] != "" : b_beta_threshold = float ( accept [ beta_key ] ) b = pars [ b_key ] cm_x = np . mean ( np . array ( x_Arai_segment ) ) cm_y = np . mean ( np . array ( y_Arai_segment ) ) a = cm_y - b * cm_x two_sigma_beta_threshold = 2 * b_beta_threshold two_sigma_slope_threshold = abs ( two_sigma_beta_threshold * b ) b1 = b + two_sigma_slope_threshold a1 = cm_y - b1 * cm_x b2 = b - two_sigma_slope_threshold a2 = cm_y - b2 * cm_x slop1 = old_div ( a1 , ( ( old_div ( a2 , b2 ) ) ) ) intercept1 = a1 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 ] ymin = intercept1 + x_Arai_segment * slop1 ymax = intercept2 + x_Arai_segment * slop2 check_1 = y_Arai_segment > ymax check_2 = y_Arai_segment < ymin if ( sum ( check_1 ) + sum ( check_2 ) ) > 0 : pars [ "fail_arai_beta_box_scatter" ] = True if len ( x_ptrm_check_for_SCAT ) > 0 : ymin = intercept1 + x_ptrm_check_for_SCAT * slop1 ymax = intercept2 + x_ptrm_check_for_SCAT * slop2 check_1 = y_ptrm_check_for_SCAT > ymax check_2 = y_ptrm_check_for_SCAT < ymin if ( sum ( check_1 ) + sum ( check_2 ) ) > 0 : pars [ "fail_ptrm_beta_box_scatter" ] = True if len ( x_tail_check_for_SCAT ) > 0 : ymin = intercept1 + x_tail_check_for_SCAT * slop1 ymax = intercept2 + x_tail_check_for_SCAT * slop2 check_1 = y_tail_check_for_SCAT > ymax check_2 = y_tail_check_for_SCAT < ymin if ( sum ( check_1 ) + sum ( check_2 ) ) > 0 : pars [ "fail_tail_beta_box_scatter" ] = True 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
calculate the paleointensity magic parameters  make some definitions