code: def eclipse_depth(mafn,Rp,Rs,b,u1=0.394,u2=0.261,max_only=False,npts=100,force_1d=False):
      
      k = Rp*REARTH/(Rs*RSUN)
      
      if max_only:
          return 1 - mafn(k,b,u1,u2)
  
      if np.size(b) == 1:
          x = np.linspace(0,np.sqrt(1-b**2),npts)
          y = b
          zs = np.sqrt(x**2 + y**2)
          fs = mafn(k,zs,u1,u2) # returns array of shape (nks,nzs)
          depth = 1-fs
      else:
          xmax = np.sqrt(1-b**2)
          x = np.linspace(0,1,npts)*xmax[:,Nones]
          y = b[:,None]
          zs = np.sqrt(x**2 + y**2)
          fs = mafn(k,zs.ravel(),u1,u2)
          if not force_1d:
              fs = fs.reshape(size(k),*zs.shape)
          depth = 1-fs
          
      meandepth = np.squeeze(depth.mean(axis=depth.ndim-1))
      
      return meandepth
docs: Calculates average (or max) eclipse depth
  
      ***why does b>1 take so freaking long?...