def York_Regression(x_segment, y_segment, x_mean, y_mean, n, lab_dc_field, steps_Arai):
    """
    input: x_segment, y_segment, x_mean, y_mean, n, lab_dc_field, steps_Arai
    output: x_err,  y_err, x_tag, y_tag, b, b_sigma, specimen_b_beta, y_intercept,
            x_intercept, x_prime, y_prime, delta_x_prime, delta_y_prime, f_Coe,
            g_Coe, g_lim, specimen_q, specimen_w, count_IZ, count_ZI, B_lab, B_anc,
            B_anc_sigma, specimen_int
    """
    x_err = x_segment - x_mean
    y_err = y_segment - y_mean
    york_b = -1* numpy.sqrt(old_div(sum(y_err**2), sum(x_err**2)) )  # averaged slope

    b = numpy.sign(sum(x_err * y_err)) * numpy.std(y_segment, ddof=1)/numpy.std(x_segment, ddof=1) # ddof is degrees of freedom
    if b == 0:
        york_b = 1e-10
    else:
        york_b = b

    york_sigma= numpy.sqrt( old_div((2 * sum(y_err**2) - 2*york_b* sum(x_err*y_err)), ( (n-2) * sum(x_err**2) )) )
    if york_sigma == 0: # prevent divide by zero
        york_sigma = 1e-10
    beta_Coe = abs(old_div(york_sigma,york_b))
    # y_T is the intercept of the extrepolated line
    # through the center of mass (see figure 7 in Coe (1978))
    y_T = y_mean - (york_b* x_mean)
    x_T = old_div((-1 * y_T), york_b)  # x intercept
    # # calculate the extrarpolated data points for f and fvds
    x_tag = old_div((y_segment - y_T ), york_b) # returns array of y points minus the y intercept, divided by slope
    y_tag = york_b*x_segment + y_T

    # intersect of the dashed square and the horizontal dahed line  next to delta-y-5 in figure 7, Coe (1978)
    x_prime = old_div((x_segment+x_tag), 2)
    y_prime = old_div((y_segment+y_tag), 2)

    delta_x_prime = abs(max(x_prime) - min(x_prime)) # TRM length of best fit line
    delta_y_prime = abs(max(y_prime) - min(y_prime)) # NRM length of best fit line

    f_Coe = old_div(delta_y_prime, abs(y_T))
    if delta_y_prime:
        g_Coe =  1 - (old_div(sum((y_prime[:-1]-y_prime[1:])**2), delta_y_prime ** 2))  # gap factor
    else:
        g_Coe = float('nan')
    g_lim = old_div((float(n) - 2), (float(n) - 1))
    q_Coe = abs(york_b)*f_Coe*g_Coe/york_sigma
    w_Coe = old_div(q_Coe, numpy.sqrt(n - 2))
    count_IZ = steps_Arai.count('IZ')
    count_ZI = steps_Arai.count('ZI')
    B_lab = lab_dc_field * 1e6
    B_anc = abs(york_b) * B_lab # in microtesla
    B_anc_sigma = york_sigma * B_lab
    specimen_int = -1* lab_dc_field * york_b # in tesla
    specimen_int_sigma = york_sigma * lab_dc_field
    return {'x_err': x_err, 'y_err': y_err, 'x_tag': x_tag, 'y_tag': y_tag,
            'specimen_b': york_b, 'specimen_b_sigma': york_sigma, 'specimen_b_beta': beta_Coe,
            'y_int': y_T, 'x_int': x_T, 'x_prime': x_prime, 'y_prime': y_prime,
            'delta_x_prime': delta_x_prime, 'delta_y_prime': delta_y_prime, 'specimen_f': f_Coe,
            'specimen_g': g_Coe, 'specimen_g_lim': g_lim, 'specimen_q': q_Coe, 'specimen_w': w_Coe,
            'count_IZ': count_IZ, 'count_ZI': count_ZI, 'B_lab': B_lab, 'B_anc': B_anc,
            'B_anc_sigma': B_anc_sigma, 'specimen_int': specimen_int, 'specimen_int_sigma': specimen_int_sigma}