Returns random eclipse parameters for provided inputs


    :param M1s, M2s, R1s, R2s, mag1s, mag2s: (array-like)
        Primary and secondary properties (mass, radius, magnitude)

    :param u11s, u21s, u12s, u22s: (optional)
        Limb darkening parameters (u11 = u1 for star 1, u21 = u2 for star 1, etc.)

    :param Ps: (array-like, optional)
        Orbital periods; same size as ``M1s``, etc.
        If only a single period is desired, use ``period``.

    :param period: (optional)
        Orbital period; use this keyword if only a single period is desired.

    :param logperkde: (optional)
        If neither ``Ps`` nor ``period`` is provided, then periods will be
        randomly generated according to this log-period distribution.
        Default is taken from the Raghavan (2010) period distribution.

    :param incs, eccs: (optional)
        Inclinations and eccentricities.  If not passed, they will be generated.
        Eccentricities will be generated according to ``ecc_fn``; inclinations
        will be randomly generated out to ``mininc``.

    :param mininc: (optional)
        Minimum inclination to generate.  Useful if you want to enhance
        efficiency by only generating mostly eclipsing, instead of mostly
        non-eclipsing systems.  If not provided and ``calc_mininc`` is
        ``True``, then this will be calculated based on inputs.

    :param calc_mininc: (optional)
        Whether to calculate ``mininc`` based on inputs.  If truly isotropic
        inclinations are desired, set this to ``False``.

    :param maxecc: (optional)
        Maximum eccentricity to generate.

    :param ecc_fn: (callable, optional)
        Orbital eccentricity generating function.  Must return ``n`` orbital
        eccentricities generated according to provided period(s)::

            eccs = ecc_fn(n,Ps)

        Defaults to :func:`stars.utils.draw_eccs`.

    :param band: (optional)
        Photometric bandpass in which eclipse is observed.

    :param return_probability_only: (optional)
        If ``True``, then will return only the average eclipse probability
        of population.

    :param return_indices: (optional)
        If ``True``, returns the indices of the original input arrays
        that the output ``DataFrame`` corresponds to.  **This behavior
        will/should be changed to just return a ``DataFrame`` of the same
        length as inputs...**

    :param MAfn: (optional)
        :class:`transit_basic.MAInterpolationFunction` object.
        If not passed, then one with default parameters will
        be created.

    :return:
        * [``wany``: indices describing which of the original input
          arrays the output ``DataFrame`` corresponds to.
        * ``df``: ``DataFrame`` with the following columns:
          ``[{band}_mag_tot, P, ecc, inc, w, dpri, dsec,
             T14_pri, T23_pri, T14_sec, T23_sec, b_pri,
             b_sec, {band}_mag_1, {band}_mag_2, fluxfrac_1,
             fluxfrac_2, switched, u1_1, u2_1, u1_2, u2_2]``.
             **N.B. that this will be shorter than your input arrays,
             because not everything will eclipse; this behavior
             will likely be changed in the future because it's confusing.**
        * ``(prob, dprob)`` Eclipse probability with Poisson uncertainty

def calculate_eclipses(M1s, M2s, R1s, R2s, mag1s, mag2s,
                       u11s=0.394, u21s=0.296, u12s=0.394, u22s=0.296,
                       Ps=None, period=None, logperkde=RAGHAVAN_LOGPERKDE,
                       incs=None, eccs=None,
                       mininc=None, calc_mininc=True,
                       maxecc=0.97, ecc_fn=draw_eccs,
                       band='Kepler',
                       return_probability_only=False, return_indices=True,
                       MAfn=None):
    """Returns random eclipse parameters for provided inputs


    :param M1s, M2s, R1s, R2s, mag1s, mag2s: (array-like)
        Primary and secondary properties (mass, radius, magnitude)

    :param u11s, u21s, u12s, u22s: (optional)
        Limb darkening parameters (u11 = u1 for star 1, u21 = u2 for star 1, etc.)

    :param Ps: (array-like, optional)
        Orbital periods; same size as ``M1s``, etc.
        If only a single period is desired, use ``period``.

    :param period: (optional)
        Orbital period; use this keyword if only a single period is desired.

    :param logperkde: (optional)
        If neither ``Ps`` nor ``period`` is provided, then periods will be
        randomly generated according to this log-period distribution.
        Default is taken from the Raghavan (2010) period distribution.

    :param incs, eccs: (optional)
        Inclinations and eccentricities.  If not passed, they will be generated.
        Eccentricities will be generated according to ``ecc_fn``; inclinations
        will be randomly generated out to ``mininc``.

    :param mininc: (optional)
        Minimum inclination to generate.  Useful if you want to enhance
        efficiency by only generating mostly eclipsing, instead of mostly
        non-eclipsing systems.  If not provided and ``calc_mininc`` is
        ``True``, then this will be calculated based on inputs.

    :param calc_mininc: (optional)
        Whether to calculate ``mininc`` based on inputs.  If truly isotropic
        inclinations are desired, set this to ``False``.

    :param maxecc: (optional)
        Maximum eccentricity to generate.

    :param ecc_fn: (callable, optional)
        Orbital eccentricity generating function.  Must return ``n`` orbital
        eccentricities generated according to provided period(s)::

            eccs = ecc_fn(n,Ps)

        Defaults to :func:`stars.utils.draw_eccs`.

    :param band: (optional)
        Photometric bandpass in which eclipse is observed.

    :param return_probability_only: (optional)
        If ``True``, then will return only the average eclipse probability
        of population.

    :param return_indices: (optional)
        If ``True``, returns the indices of the original input arrays
        that the output ``DataFrame`` corresponds to.  **This behavior
        will/should be changed to just return a ``DataFrame`` of the same
        length as inputs...**

    :param MAfn: (optional)
        :class:`transit_basic.MAInterpolationFunction` object.
        If not passed, then one with default parameters will
        be created.

    :return:
        * [``wany``: indices describing which of the original input
          arrays the output ``DataFrame`` corresponds to.
        * ``df``: ``DataFrame`` with the following columns:
          ``[{band}_mag_tot, P, ecc, inc, w, dpri, dsec,
             T14_pri, T23_pri, T14_sec, T23_sec, b_pri,
             b_sec, {band}_mag_1, {band}_mag_2, fluxfrac_1,
             fluxfrac_2, switched, u1_1, u2_1, u1_2, u2_2]``.
             **N.B. that this will be shorter than your input arrays,
             because not everything will eclipse; this behavior
             will likely be changed in the future because it's confusing.**
        * ``(prob, dprob)`` Eclipse probability with Poisson uncertainty

    """
    if MAfn is None:
        logging.warning('MAInterpolationFunction not passed, so generating one...')
        MAfn = MAInterpolationFunction(nzs=200,nps=400,pmin=0.007,pmax=1/0.007)

    M1s = np.atleast_1d(M1s)
    M2s = np.atleast_1d(M2s)
    R1s = np.atleast_1d(R1s)
    R2s = np.atleast_1d(R2s)

    nbad = (np.isnan(M1s) | np.isnan(M2s) | np.isnan(R1s) | np.isnan(R2s)).sum()
    if nbad > 0:
        logging.warning('{} M1s are nan'.format(np.isnan(M1s).sum()))
        logging.warning('{} M2s are nan'.format(np.isnan(M2s).sum()))
        logging.warning('{} R1s are nan'.format(np.isnan(R1s).sum()))
        logging.warning('{} R2s are nan'.format(np.isnan(R2s).sum()))

    mag1s = mag1s * np.ones_like(M1s)
    mag2s = mag2s * np.ones_like(M1s)
    u11s = u11s * np.ones_like(M1s)
    u21s = u21s * np.ones_like(M1s)
    u12s = u12s * np.ones_like(M1s)
    u22s = u22s * np.ones_like(M1s)

    n = np.size(M1s)

    #a bit clunky here, but works.
    simPs = False
    if period:
        Ps = np.ones(n)*period
    else:
        if Ps is None:
            Ps = 10**(logperkde.rvs(n))
            simPs = True
    simeccs = False
    if eccs is None:
        if not simPs and period is not None:
            eccs = ecc_fn(n,period,maxecc=maxecc)
        else:
            eccs = ecc_fn(n,Ps,maxecc=maxecc)
        simeccs = True

    bad_Ps = np.isnan(Ps)
    if bad_Ps.sum()>0:
        logging.warning('{} nan periods.  why?'.format(bad_Ps.sum()))
    bad_eccs = np.isnan(eccs)
    if bad_eccs.sum()>0:
        logging.warning('{} nan eccentricities.  why?'.format(bad_eccs.sum()))

    semimajors = semimajor(Ps, M1s+M2s)*AU #in AU

    #check to see if there are simulated instances that are
    # too close; i.e. periastron sends secondary within roche
    # lobe of primary
    tooclose = withinroche(semimajors*(1-eccs)/AU,M1s,R1s,M2s,R2s)
    ntooclose = tooclose.sum()
    tries = 0
    maxtries=5
    if simPs:
        while ntooclose > 0:
            lastntooclose=ntooclose
            Ps[tooclose] = 10**(logperkde.rvs(ntooclose))
            if simeccs:
                eccs[tooclose] = draw_eccs(ntooclose,Ps[tooclose])
            semimajors[tooclose] = semimajor(Ps[tooclose],M1s[tooclose]+M2s[tooclose])*AU
            tooclose = withinroche(semimajors*(1-eccs)/AU,M1s,R1s,M2s,R2s)
            ntooclose = tooclose.sum()
            if ntooclose==lastntooclose:   #prevent infinite loop
                tries += 1
                if tries > maxtries:
                    logging.info('{} binaries are "too close"; gave up trying to fix.'.format(ntooclose))
                    break
    else:
        while ntooclose > 0:
            lastntooclose=ntooclose
            if simeccs:
                eccs[tooclose] = draw_eccs(ntooclose,Ps[tooclose])
            semimajors[tooclose] = semimajor(Ps[tooclose],M1s[tooclose]+M2s[tooclose])*AU
            #wtooclose = where(semimajors*(1-eccs) < 2*(R1s+R2s)*RSUN)
            tooclose = withinroche(semimajors*(1-eccs)/AU,M1s,R1s,M2s,R2s)
            ntooclose = tooclose.sum()
            if ntooclose==lastntooclose:   #prevent infinite loop
                tries += 1
                if tries > maxtries:
                    logging.info('{} binaries are "too close"; gave up trying to fix.'.format(ntooclose))
                    break

    #randomize inclinations, either full range, or within restricted range
    if mininc is None and calc_mininc:
        mininc = minimum_inclination(Ps, M1s, M2s, R1s, R2s)

    if incs is None:
        if mininc is None:
            incs = np.arccos(np.random.random(n)) #random inclinations in radians
        else:
            incs = np.arccos(np.random.random(n)*np.cos(mininc*np.pi/180))
    if mininc:
        prob = np.cos(mininc*np.pi/180)
    else:
        prob = 1

    logging.debug('initial probability given mininc starting at {}'.format(prob))

    ws = np.random.random(n)*2*np.pi

    switched = (R2s > R1s)
    R_large = switched*R2s + ~switched*R1s
    R_small = switched*R1s + ~switched*R2s


    b_tras = semimajors*np.cos(incs)/(R_large*RSUN) * (1-eccs**2)/(1 + eccs*np.sin(ws))
    b_occs = semimajors*np.cos(incs)/(R_large*RSUN) * (1-eccs**2)/(1 - eccs*np.sin(ws))

    b_tras[tooclose] = np.inf
    b_occs[tooclose] = np.inf

    ks = R_small/R_large
    Rtots = (R_small + R_large)/R_large
    tra = (b_tras < Rtots)
    occ = (b_occs < Rtots)
    nany = (tra | occ).sum()
    peb = nany/float(n)
    prob *= peb
    if return_probability_only:
        return prob,prob*np.sqrt(nany)/n


    i = (tra | occ)
    wany = np.where(i)
    P,M1,M2,R1,R2,mag1,mag2,inc,ecc,w = Ps[i],M1s[i],M2s[i],R1s[i],R2s[i],\
        mag1s[i],mag2s[i],incs[i]*180/np.pi,eccs[i],ws[i]*180/np.pi
    a = semimajors[i]  #in cm already
    b_tra = b_tras[i]
    b_occ = b_occs[i]
    u11 = u11s[i]
    u21 = u21s[i]
    u12 = u12s[i]
    u22 = u22s[i]


    switched = (R2 > R1)
    R_large = switched*R2 + ~switched*R1
    R_small = switched*R1 + ~switched*R2
    k = R_small/R_large

    #calculate durations
    T14_tra = P/np.pi*np.arcsin(R_large*RSUN/a * np.sqrt((1+k)**2 - b_tra**2)/np.sin(inc*np.pi/180)) *\
        np.sqrt(1-ecc**2)/(1+ecc*np.sin(w*np.pi/180)) #*24*60
    T23_tra = P/np.pi*np.arcsin(R_large*RSUN/a * np.sqrt((1-k)**2 - b_tra**2)/np.sin(inc*np.pi/180)) *\
        np.sqrt(1-ecc**2)/(1+ecc*np.sin(w*np.pi/180)) #*24*60
    T14_occ = P/np.pi*np.arcsin(R_large*RSUN/a * np.sqrt((1+k)**2 - b_occ**2)/np.sin(inc*np.pi/180)) *\
        np.sqrt(1-ecc**2)/(1-ecc*np.sin(w*np.pi/180)) #*24*60
    T23_occ = P/np.pi*np.arcsin(R_large*RSUN/a * np.sqrt((1-k)**2 - b_occ**2)/np.sin(inc*np.pi/180)) *\
        np.sqrt(1-ecc**2)/(1-ecc*np.sin(w*np.pi/180)) #*24*60

    bad = (np.isnan(T14_tra) & np.isnan(T14_occ))
    if bad.sum() > 0:
        logging.error('Something snuck through with no eclipses!')
        logging.error('k: {}'.format(k[bad]))
        logging.error('b_tra: {}'.format(b_tra[bad]))
        logging.error('b_occ: {}'.format(b_occ[bad]))
        logging.error('T14_tra: {}'.format(T14_tra[bad]))
        logging.error('T14_occ: {}'.format(T14_occ[bad]))
        logging.error('under sqrt (tra): {}'.format((1+k[bad])**2 - b_tra[bad]**2))
        logging.error('under sqrt (occ): {}'.format((1+k[bad])**2 - b_occ[bad]**2))
        logging.error('eccsq: {}'.format(ecc[bad]**2))
        logging.error('a in Rsun: {}'.format(a[bad]/RSUN))
        logging.error('R_large: {}'.format(R_large[bad]))
        logging.error('R_small: {}'.format(R_small[bad]))
        logging.error('P: {}'.format(P[bad]))
        logging.error('total M: {}'.format(M1[bad]+M2[bad]))

    T14_tra[(np.isnan(T14_tra))] = 0
    T23_tra[(np.isnan(T23_tra))] = 0
    T14_occ[(np.isnan(T14_occ))] = 0
    T23_occ[(np.isnan(T23_occ))] = 0

    #calling mandel-agol
    ftra = MAfn(k,b_tra,u11,u21)
    focc = MAfn(1/k,b_occ/k,u12,u22)

    #fix those with k or 1/k out of range of MAFN....or do it in MAfn eventually?
    wtrabad = np.where((k < MAfn.pmin) | (k > MAfn.pmax))
    woccbad = np.where((1/k < MAfn.pmin) | (1/k > MAfn.pmax))
    for ind in wtrabad[0]:
        ftra[ind] = occultquad(b_tra[ind],u11[ind],u21[ind],k[ind])
    for ind in woccbad[0]:
        focc[ind] = occultquad(b_occ[ind]/k[ind],u12[ind],u22[ind],1/k[ind])

    F1 = 10**(-0.4*mag1) + switched*10**(-0.4*mag2)
    F2 = 10**(-0.4*mag2) + switched*10**(-0.4*mag1)

    dtra = 1-(F2 + F1*ftra)/(F1+F2)
    docc = 1-(F1 + F2*focc)/(F1+F2)

    totmag = -2.5*np.log10(F1+F2)

    #wswitched = where(switched)
    dtra[switched],docc[switched] = (docc[switched],dtra[switched])
    T14_tra[switched],T14_occ[switched] = (T14_occ[switched],T14_tra[switched])
    T23_tra[switched],T23_occ[switched] = (T23_occ[switched],T23_tra[switched])
    b_tra[switched],b_occ[switched] = (b_occ[switched],b_tra[switched])
    #mag1[wswitched],mag2[wswitched] = (mag2[wswitched],mag1[wswitched])
    F1[switched],F2[switched] = (F2[switched],F1[switched])
    u11[switched],u12[switched] = (u12[switched],u11[switched])
    u21[switched],u22[switched] = (u22[switched],u21[switched])

    dtra[(np.isnan(dtra))] = 0
    docc[(np.isnan(docc))] = 0

    if np.any(np.isnan(ecc)):
        logging.warning('{} nans in eccentricity.  why?'.format(np.isnan(ecc).sum()))

    df =  pd.DataFrame({'{}_mag_tot'.format(band) : totmag,
                        'P':P, 'ecc':ecc, 'inc':inc, 'w':w,
                        'dpri':dtra, 'dsec':docc,
                        'T14_pri':T14_tra, 'T23_pri':T23_tra,
                        'T14_sec':T14_occ, 'T23_sec':T23_occ,
                        'b_pri':b_tra, 'b_sec':b_occ,
                        '{}_mag_1'.format(band) : mag1,
                        '{}_mag_2'.format(band) : mag2,
                        'fluxfrac_1':F1/(F1+F2),
                        'fluxfrac_2':F2/(F1+F2),
                        'switched':switched,
                        'u1_1':u11, 'u2_1':u21, 'u1_2':u12, 'u2_2':u22})

    df.reset_index(inplace=True)

    logging.debug('final prob: {}'.format(prob))

    if return_indices:
        return wany, df, (prob, prob*np.sqrt(nany)/n)
    else:
        return df, (prob, prob*np.sqrt(nany)/n)