forked from SgmAstro/DESI
-
Notifications
You must be signed in to change notification settings - Fork 1
/
ddp.py
93 lines (60 loc) · 3.28 KB
/
ddp.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
import os
import fitsio
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import Table
from cosmo import volcom
from scipy.interpolate import interp1d
from findfile import findfile
tmr_DDP1 = [-21.8, -20.1]
tmr_DDP2 = [-20.6, -19.3]
tmr_DDP3 = [-19.6, -17.8]
def initialise_ddplimits(survey, Mcol='M0P0_QALL', bright_idx=None, faint_idx=None):
if bright_idx == None:
bright_idx = 3
assert Mcol == 'M0P0_QALL'
if faint_idx == None:
faint_idx = 17
assert Mcol == 'M0P0_QALL'
return _initialise_ddplimits(bright_idx, faint_idx, survey=survey)
def _initialise_ddplimits(bright_idx, faint_idx, Mcol='M0P0_QALL', survey='gama'):
bpath = findfile(ftype='ddp_limit', dryrun=False, survey=survey, ddp_count=bright_idx)
fpath = findfile(ftype='ddp_limit', dryrun=False, survey=survey, ddp_count=faint_idx)
print(f'Reading {bpath}')
print(f'Reading {fpath}')
_bright_curve = Table.read(bpath)
_faint_curve = Table.read(fpath)
# TODO: extend the curve limits and put bounds_error back on.
bright_curve = interp1d(_bright_curve[Mcol], _bright_curve['Z'], kind='linear', copy=True, bounds_error=False, fill_value=0.0, assume_sorted=False)
bright_curve_r = interp1d(_bright_curve['Z'], _bright_curve[Mcol], kind='linear', copy=True, bounds_error=False, fill_value=0.0, assume_sorted=False)
faint_curve = interp1d(_faint_curve[Mcol], _faint_curve['Z'], kind='linear', copy=True, bounds_error=False, fill_value=1.0, assume_sorted=False)
faint_curve_r = interp1d(_faint_curve['Z'], _faint_curve[Mcol], kind='linear', copy=True, bounds_error=False, fill_value=1.0, assume_sorted=False)
return bright_curve, bright_curve_r, faint_curve, faint_curve_r
def get_ddps(Area, M_0P0s, zs, survey):
result = np.zeros(len(zs) * 3, dtype=int).reshape(len(zs), 3)
resultz = np.zeros(len(zs) * 3, dtype=int).reshape(len(zs), 3)
zlims = {}
bright_curve, bright_curve_r, faint_curve, faint_curve_r = initialise_ddplimits(survey=survey)
for i, lims in enumerate([tmr_DDP1, tmr_DDP2, tmr_DDP3]):
in_ddp = (M_0P0s >= lims[0]) & (M_0P0s <= lims[1])
zmax = np.atleast_1d(faint_curve(lims[1]))[0]
zmin = np.atleast_1d(bright_curve(lims[0]))[0]
exclude = (zs > zmax) | (zs < zmin)
in_ddp = in_ddp & ~exclude
in_ddpz = ~exclude
result[in_ddp, i] = 1
resultz[in_ddpz, i] = 1
ddp_zs = zs[in_ddp]
print(zmin, zmax, len(ddp_zs))
zmax = np.array([zmax, ddp_zs.max()]).min()
zmin = np.array([zmin, ddp_zs.min()]).max()
zlims['DDP{}_ZMIN'.format(i+1)] = zmin
zlims['DDP{}_ZMAX'.format(i+1)] = zmax
zlims['DDP{}_VZ'.format(i+1)] = volcom(zmax, Area) - volcom(zmin, Area)
zlims['DDP{}ZLIMS_NGAL'.format(i+1)] = np.count_nonzero(in_ddpz)
zlims['DDP{}_NGAL'.format(i+1)] = np.count_nonzero(in_ddp)
zlims['DDP{}_DENS'.format(i+1)] = np.count_nonzero(in_ddp) / zlims['DDP{}_VZ'.format(i+1)]
return result, resultz, zlims
if __name__ == '__main__':
initialise_ddplimits('desi', Mcol='M0P0_QALL')
print('Done.')