add first version of code

FFD.py

@@ -0,0 +1,140 @@
+import numpy as np
+from astropy.modeling.models import Gaussian2D
+
+
+def FFD(ED, TOTEXP=1., Lum=30., fluxerr=0., dur=[], logY=True, est_comp=False):
+    '''
+    Given a set of stellar flares, with accompanying durations light curve properties,
+    compute the reverse cumulative Flare Frequency Distribution (FFD), and
+    approximate uncertainties in both energy and rate (X,Y) dimensions.
+    This diagram can be read as measuring the number of flares per day at a
+    given energy or larger.
+
+    Not a complicated task, just tedious.
+
+    Y-errors (rate) are computed using Poisson upper-limit approximation from
+    Gehrels (1986) "Confidence limits for small numbers of events in astrophysical data", https://doi.org/10.1086/164079
+    Eqn 7, assuming S=1.
+
+    X-errors (event energy) are computed following Signal-to-Noise approach commonly
+    used for Equivalent Widths in spectroscopy, from
+    Vollmann & Eversberg (2006) "Astronomische Nachrichten, Vol.327, Issue 9, p.862", https://dx.doi.org/10.1002/asna.200610645
+    Eqn 6.
+
+    Parameters
+    ----------
+    ED : array of Equiv Dur's, need to include a luminosity!
+    TOTEXP : total duration of observations, in days
+    Lum : the luminosity of the star
+    fluxerr : the average flux errors of your data (in relative flux units!)
+    dur : array of flare durations.
+    logY : if True return Y-axis (and error) in log rate (Default: True)
+    est_comp : estimate incompleteness using histogram method, scale Y errors?
+        (Default: True)
+
+    Returns
+    -------
+    ffd_x, ffd_y, ffd_xerr, ffd_yerr
+
+    X coordinate always assumed to be log_10(Energy)
+    Y coordinate is log_10(N/Day) by default, but optionally is N/Day
+
+    Upgrade Ideas
+    -------------
+    - More graceful behavior if only an array of flares and a total duration are
+        specified (i.e. just enough to make ffd_x, ffd_y)
+    - Better propogation of specific flux errors in the light curve, rather than
+        average error used
+    - Include detrending errors? (e.g. from a GP)
+    - Asymmetric Poisson errors?
+    - Better handling of incompleteness?
+
+    '''
+    # REVERSE sort the flares in energy
+    ss = np.argsort(np.array(ED))[::-1]
+    ffd_x = np.log10(ED[ss]) + Lum
+
+    Num = np.arange(1, len(ffd_x)+1)
+    ffd_y = Num / TOTEXP
+
+    # approximate the Poisson Y errors using Gehrels (1986) eqn 7
+    Perror = np.sqrt(Num + 0.75) + 1.0
+    ffd_yerr = Perror / TOTEXP
+
+    # estimate completeness using the cumulative distribution of the histogram
+    if est_comp:
+        # make very loose guess at how many bins to choose
+        nbin = int(np.sqrt(len(ffd_x)))
+        if nbin < 10:
+            nbin=10 # but use at least 10 bins
+
+        # make histogram of the log(energies)
+        hh, be = np.histogram(ffd_x, bins=nbin, range=[np.nanmin(ffd_x), np.nanmax(ffd_x)])
+        hh = hh/np.nanmax(hh)
+        # make cumulative distribution of the histogram, scale to =1 at the hist peak
+        cc = np.cumsum(hh)/np.sum(hh[0:np.argmax(hh)])
+        be = (be[1:]+be[0:-1])/2
+        # make completeness = 1 for energies above the histogram peak
+        cc[np.argmax(hh):] = 1
+        # interpolate the cumulative histogram curve back to the original energies
+        ycomp = np.interp(ffd_x, be, cc)
+        # scale the y-errors by the completeness factor (i.e. inflate small energy errors)
+        ffd_yerr = ffd_yerr / ycomp
+
+    if logY:
+        # transform FFD Y and Y Error into log10
+        ffd_yerr = np.abs(ffd_yerr / np.log(10.) / ffd_y)
+        ffd_y = np.log10(ffd_y)
+
+    # compute X uncertainties for FFD
+    if len(dur)==len(ffd_x):
+
+        # assume relative flux error = 1/SN
+        S2N = 1/fluxerr
+        # based on Equivalent Width error
+        # Eqn 6, Vollmann & Eversberg (2006) Astronomische Nachrichten, Vol.327, Issue 9, p.862
+        ED_err = np.sqrt(2)*(dur[ss]*86400. - ED[ss])/S2N
+        ffd_xerr = np.abs((ED_err) / np.log(10.) / ED[ss]) # convert to log
+    else:
+        # not particularly meaningful, but an interesting shape. NOT reccomended
+        print('Warning: Durations not set. Making bad assumptions about the FFD X Error!')
+        ffd_xerr = (1/np.sqrt(ffd_x-np.nanmin(xT))/(np.nanmax(ffd_x)-np.nanmin(ffd_x)))
+
+    return ffd_x, ffd_y, ffd_xerr, ffd_yerr
+
+
+def FlareKernel(x, y, xe, ye, Nx=100, Ny=100, xlim=[], ylim=[], return_axis=True):
+    '''
+    Use 2D Gaussians (from astropy models) to make a basic kernel density,
+    with errors in both X and Y considered. Turn into a 2D "image"
+
+    Upgrade Ideas
+    -------------
+    It's slow. Since Gaussians are defined analytically, maybe this could be
+    re-cast as a single array math opperation, and then refactored to have the
+    same fit/evaluate behavior as KDE functions.  Hmm...
+    '''
+
+    if len(xlim) == 0:
+        xlim = [np.nanmin(x) - np.nanmean(xe), np.nanmax(x) + np.nanmean(xe)]
+    if len(ylim) == 0:
+        ylim = [np.nanmin(y) - np.nanmean(ye), np.nanmax(y) + np.nanmean(ye)]
+
+    xx,yy = np.meshgrid(np.linspace(xlim[0], xlim[1], Nx),
+                        np.linspace(ylim[0], ylim[1], Ny), indexing='xy')
+    dx = (np.max(xlim)-np.min(xlim)) / (Nx-1)
+    dy = (np.max(ylim)-np.min(ylim)) / (Ny-1)
+
+    im = np.zeros_like(xx)
+
+    for k in range(len(x)):
+        g = Gaussian2D(amplitude=1/(2*np.pi*(xe[k]+dx)*(ye[k]+dy)),
+                       x_mean=x[k], y_mean=y[k], x_stddev=xe[k]+dx, y_stddev=ye[k]+dy)
+        tmp = g(xx,yy)
+        if np.isfinite(np.sum(tmp)):
+            im = im + tmp
+
+    if return_axis:
+        return im, xx, yy
+    else:
+        return im

README.md

@@ -1,2 +1,4 @@
 # FFD
 How to generate a stellar Flare Frequency Distribution, and its uncertainties
+
+<img src="https://github.com/jradavenport/FFD/blob/master/ffd.png" alt="ffd package example" width="400"/>

__init__.py

@@ -0,0 +1,4 @@
+__version__ = "0.1"
+__author__ = "James Davenport (@jradavenport)"
+
+from .FFD import (FFD, FlareKernel)