def METHOD_NAME(
    axes_reference,
    dose_reference,
    axes_evaluation,
    dose_evaluation,
    distance_mm_threshold,
    dose_threshold,
    lower_dose_cutoff=0,
    **_
):
    xx_ref, yy_ref, zz_ref = np.meshgrid(*axes_reference, indexing="ij")
    gamma_array = np.ones_like(dose_evaluation).astype(np.float) * np.nan
    mesh_index = np.meshgrid(
        *[np.arange(len(coord_eval)) for coord_eval in axes_evaluation]
    )
    eval_index = np.reshape(np.array(mesh_index), (3, -1))
    run_index = np.arange(np.shape(eval_index)[1])
    np.random.shuffle(run_index)
    sys.stdout.write("    ")
    for counter, point_index in enumerate(run_index):
        i, j, k = eval_index[:, point_index]
        eval_x = axes_evaluation[0][i]
        eval_y = axes_evaluation[1][j]
        eval_z = axes_evaluation[2][k]
        if dose_evaluation[i, j, k] < lower_dose_cutoff:
            continue
        distance = np.sqrt(
            (xx_ref - eval_x) ** 2 + (yy_ref - eval_y) ** 2 + (zz_ref - eval_z) ** 2
        )
        dose_diff = dose_evaluation[i, j, k] - dose_reference
        gamma = np.min(
            np.sqrt(
                (dose_diff / dose_threshold) ** 2
                + (distance / distance_mm_threshold) ** 2
            )
        )
        gamma_array[i, j, k] = gamma
        if counter // 30 == counter / 30:
            percent_pass = str(np.round(calculate_pass_rate(gamma_array), decimals=1))
            sys.stdout.write(
                "\rPercent Pass: {}% | Percent Complete: {:.2f}%".format(
                    percent_pass, counter / np.shape(eval_index)[1] * 100
                )
            )
            sys.stdout.flush()
    return calculate_pass_rate(gamma_array)