Number_Counts/Number_Counts.py at master · JoshGreenslade/Number_Counts · GitHub

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
141
142
143
144
145
146
# -*- coding: utf-8 -*-
"""
Created on Fri Oct  9 12:40:22 2015

@author: JG1714

This is a section of code that will estimate the differential and integral
number counts of sources from a catalogue.

In general the catalogue should only need the list of fluxes.

Before passing the list,
"""

# =============================================================================
# Integral counts of the blank field number sources (schecter function)

nq = 1599
sp = 3.3
alpha = -2.0

# Define the schecter function
schec = lambda s: (nq / sp) * s * (np.power((s / sp),alpha)) * np.exp(-s / sp)

fluxes = np.array([])
n = 1
ma = 30
mi = 1

integral_counts = integrate.quad(schec, mi, ma)
s_max = schec(mi)
s_min = schec(ma)
while np.size(fluxes) < int(integral_counts[0]*0.2):
    n += 1
    test =  np.random.rand()*np.abs(s_max-s_min) + s_min
    test2 = np.random.rand()*np.abs(ma-mi)+mi


    if test <= schec(test2):
        fluxes = np.append(fluxes, test2)
        print(n)


fluxes = fluxes #### !!! Replace this line with the flux values !!! ###
area = 0.2 ### !!! Replace this line with the area (in square degrees) of survey.


####################### Integral Number Counts #############################

"""
This section will calculate the integral number counts from the catalogue.
"""

# Input parameters to the system

flux_minimum = 2   # The minimum flux to count from
flux_maximum = 20    # The maximum flux to count up to
n_bins = 30          # The number of bins to use


"""
For each bin, calculate the number of sources in it
"""

# Set up the bins
bin_edges = np.linspace(flux_minimum,
                        flux_maximum,
                        n_bins + 1)  # As we will end up with n-1 bin centres
bin_centres = ( bin_edges + np.roll(bin_edges, -1) )/2
bin_centres = bin_centres[:-1]
bin_counts = np.zeros(n_bins) # n_bins - 1 as only values between values


# For each bin, add all the counts of sources within that bin

for index, Bin in enumerate(bin_edges[:-1]):

    n_sources = np.sum(fluxes >= Bin)
    bin_counts[index] = n_sources


# The errors on each bin are poisson in nature, so sqrt for error

bin_errors = np.sqrt(bin_counts)


# Convert these to the raw bin counts, so we can modify them later yet retain
# their values.

raw_Ibin_counts = bin_counts
raw_Ibin_errors = bin_errors

raw_int_counts = bin_counts/area
raw_int_errors = bin_errors/area
raw_int_centres = bin_centres

"""
We now have a simple count of the number of sources that have fluxes greater
than each value in bin_edges, and an estimate of the poisson error on each of
them.

"""


####################### Differential Number Counts #############################

"""
This section will calculate the differential number counts from the catalogue.
Very similar to the previous section
"""


bin_widths = (np.roll(bin_edges, -1) - bin_edges)
bin_widths = bin_widths[:-1]

bin_counts = np.zeros(n_bins) # n_bins - 1 as only values between values


# For each bin, add all the counts of sources within that bin

for index, Bin in enumerate(bin_edges[:-1]):

    source_selection = (fluxes >= bin_edges[index]) * (fluxes < bin_edges[index+1])
    n_sources = np.sum(source_selection)
    bin_counts[index] = n_sources


# The errors on each bin are poisson in nature, so sqrt for error

bin_errors = np.sqrt(bin_counts)


# Convert these to the raw bin counts, so we can modify them later yet retain
# their values.

raw_Dbin_counts = bin_counts/bin_widths
raw_Dbin_errors = bin_errors/bin_widths

raw_dif_centres = bin_centres
raw_dif_counts = raw_Dbin_counts/area
raw_dif_errors = raw_Dbin_errors/area