Returns m1, m2, and m3 appropriate for TripleStarPopulation, given "primary" mass (most massive of system) and binary/triple fractions.


    star with m1 orbits (m2 + m3).  This means that the primary mass mA will correspond
    either to m1 or m2.  Any mass set to 0 means that component does not exist.

def mult_masses(mA, f_binary=0.4, f_triple=0.12,
                minmass=0.11, qmin=0.1, n=1e5):
    """Returns m1, m2, and m3 appropriate for TripleStarPopulation, given "primary" mass (most massive of system) and binary/triple fractions.


    star with m1 orbits (m2 + m3).  This means that the primary mass mA will correspond
    either to m1 or m2.  Any mass set to 0 means that component does not exist.
    """

    if np.size(mA) > 1:
        n = len(mA)
    else:
        mA = np.ones(n) * mA


    r = rand.random(n)
    is_single = r > (f_binary + f_triple)
    is_double = (r > f_triple) & (r < (f_binary + f_triple))
    is_triple = r <= f_triple

    CwA = rand.random(n) < 0.5
    CwB = ~CwA


    #these for Triples:

    minq2_A = minmass/mA
    q2_A = rand.random(n)*(1-minq2_A) + minq2_A
    minq1_A = (minmass/mA)/(1+q2_A)
    maxq1_A = 1/(1+q2_A)
    q1_A = rand.random(n)*(maxq1_A-minq1_A) + minq1_A

    minq1_B = 2*minmass/mA
    q1_B = rand.random(n)*(1-minq1_B) + minq1_B
    minq2_B = np.maximum(((q1_B*mA)-minmass)/minmass,
                         (q1_B*mA - minmass)/(q1_B*mA + minmass))
    maxq2_B = 1.
    q2_B = rand.random(n)*(maxq2_B-minq2_B) + minq2_B


    mB_A = q1_A*(1 + q2_A) * mA
    mC_A = q2_A * mA

    mB_B = (q1_B/(1 + q2_B)) * mA
    mC_B = (q1_B*q2_B)/(1 + q2_B) * mA

    mB = CwA*mB_A + CwB*mB_B
    mC = CwA*mC_A + CwB*mC_B

    #for binaries-only
    qmin = minmass/mA
    q = rand.random(n)*(1-qmin) + qmin
    mB[is_double] = q[is_double]*mA[is_double]


    #now need to define the proper mapping from A,B,C to 1,2,3:
    # If no B or C present, then A=1
    # If B present but not C, then A=1, B=2
    # If both B and C present then:
    #     If C is with A, then A=2, C=3, B=1
    #     If C is with B, then A=1, B=2, C=3
    m1 = (mA)*(is_single + is_double) + (CwA*mB + CwB*mA)*is_triple
    m2 = (mB)*is_double + (CwA*mA + CwB*mB)*is_triple
    m3 = mC*is_triple

    return m1, m2, m3