docs: 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)