-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathFFD.py
140 lines (114 loc) · 5.48 KB
/
FFD.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
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 Equivalent Durations of flares (units of seconds)
TOTEXP : total duration of observations (units of days)
Lum : log luminosity of the star (units of log10[erg/s])
fluxerr : the average flux errors of your data (in relative flux units)
dur : array of flare durations (units of days)
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