def dt_pi_constraint_preserving_spherical_radiation(face_mesh_velocity, normal_covector, normal_vector, psi, phi, inertial_coords, gamma1, gamma2, lapse, shift, dt_psi, dt_pi, dt_phi, d_psi, d_phi):
    dt_psi_correction = dt_psi_constraint_preserving_spherical_radiation(face_mesh_velocity, normal_covector, normal_vector, psi, phi, inertial_coords, gamma1, gamma2, lapse, shift, dt_psi, dt_pi, dt_phi, d_psi, d_phi)
    inv_radius = (1.0 / np.linalg.norm(inertial_coords))
    bc_dt_pi = (((((((2.0 * (inv_radius ** 2)) * psi) + ((4.0 * inv_radius) * dt_psi)) + ((4.0 * inv_radius) * np.dot(normal_vector, phi))) + (2.0 * np.dot(normal_vector, dt_phi))) + np.dot(shift, dt_phi)) + np.einsum('i,j,ij', normal_vector, normal_vector, d_phi))
    bc_dt_pi /= lapse
    return ((bc_dt_pi - dt_pi) + (gamma2 * dt_psi_correction))
def dt_pi_constraint_preserving_spherical_radiation(face_mesh_velocity, normal_covector, normal_vector, psi, phi, inertial_coords, gamma1, gamma2, lapse, shift, dt_psi, dt_pi, dt_phi, d_psi, d_phi):
    dt_psi_correction = dt_psi_constraint_preserving_spherical_radiation(face_mesh_velocity, normal_covector, normal_vector, psi, phi, inertial_coords, gamma1, gamma2, lapse, shift, dt_psi, dt_pi, dt_phi, d_psi, d_phi)
    inv_radius = (1.0 / np.linalg.norm(inertial_coords))
    bc_dt_pi = (((((((2.0 * (inv_radius ** 2)) * psi) + ((4.0 * inv_radius) * dt_psi)) + ((4.0 * inv_radius) * np.dot(normal_vector, phi))) + (2.0 * np.dot(normal_vector, dt_phi))) + np.dot(shift, dt_phi)) + np.einsum('i,j,ij', normal_vector, normal_vector, d_phi))
    bc_dt_pi /= lapse
    return ((bc_dt_pi - dt_pi) + (gamma2 * dt_psi_correction))