def EMD_2v3_body_fn(eta, xi, s, t):
    # Getting the angle between quark and anti-quark (qbar), and quark and photon (y) from kinematics, following mathematica calculations
    theta_qbar = np.arccos( ( 1.-s-t-s*t )/( (1.-s)*(1.-t) ) )
    theta_y    = np.arcsin( -1.*(1.-t)*np.sin(theta_qbar(s,t)) / (s+t) )

    # Defining the thrust w.r.t. a thrust axis, minimizing, and finding the thrust axis:
    thrust = lambda phi: (1.-s)*abs( np.cos(phi) ) + (1.-t)*abs( np.cos(theta_qbar(s,t)-phi) ) + (s+t)*abs( np.cos(theta_y(s,t)-phi) )
    theta_th = minimize(thrust, x0=0., bounds=((0.,np.pi),) ).x

    # Defining an EMD with 2 free parameters under which we will minimize
    EMD_fn = (
              (eta+xi-s)*np.sqrt(1.-np.cos(theta_th)) + (1.-eta-xi)*np.sqrt(1.+np.cos(theta_th))
              (1.-t-eta)*np.sqrt(1.-np.cos(theta_th-theta_qbar(s,t))) + eta*np.sqrt(1.+np.cos(theta_th-theta_qbar(s,t)))
              (s+t-xi)*np.sqrt(1.-np.cos(theta_th-theta_y(s,t))) + xi*np.sqrt(1.+np.cos(theta_th-theta_y(s,t)))
             )/np.sqrt(2)

    return EMD_fn



def MC_thrust_aligned_EMD_Distribution_Plot(numpoints):
    numSamplePoints = 100000
    # 100000 MC sample points
    st, area = sample_stSpace(1000000)

    rho_EMD = []
    error = []
    EMD_list = []

    maxEMD = 1
    EMD_vals = np.linspace(0, maxEMD, numpoints)
    epsilon = maxEMD/numpoints


    # Constraints on the solutions space, requiring that the EMD 'matrix elements' are positive
    EMD_cons = [
                {'type': 'ineq', 'fun': lambda data: eta+xi-s},
                {'type': 'ineq', 'fun': lambda data: 1.-eta-xi},
                {'type': 'ineq', 'fun': lambda data: 1.-t-eta},
                {'type': 'ineq', 'fun': lambda data: eta},
                {'type': 'ineq', 'fun': lambda data: s+t-xi},
                {'type': 'ineq', 'fun': lambda data: xi}
               ]

    # Initial guess and bounds for eta and xi
    initial_guess = (.5,.5)
    bnds = ((0,1),(0,1))

    print(minimize(EMD_fn, x0=initial_guess, bounds=bnds, constraints=cons))