1 /*
2 * Copyright (c) 2012-2021, The OSKAR Developers.
3 * See the LICENSE file at the top-level directory of this distribution.
4 */
5
6 #include "sky/oskar_sky.h"
7 #include "utility/oskar_device.h"
8
9 #include <float.h>
10
11 #ifdef __cplusplus
12 extern "C" {
13 #endif
14
oskar_sky_filter_by_flux(oskar_Sky * sky,double min_I,double max_I,int * status)15 void oskar_sky_filter_by_flux(oskar_Sky* sky, double min_I, double max_I,
16 int* status)
17 {
18 int in = 0, out = 0;
19 if (*status) return;
20
21 /* Return immediately if no filtering should be done. */
22 if (min_I <= -DBL_MAX && max_I >= DBL_MAX) return;
23
24 if (max_I < min_I)
25 {
26 *status = OSKAR_ERR_INVALID_ARGUMENT;
27 return;
28 }
29
30 /* Get the meta-data. */
31 const int location = oskar_sky_mem_location(sky);
32 const int type = oskar_sky_precision(sky);
33 const int num_sources = oskar_sky_num_sources(sky);
34
35 /* Filtering is only supported for data in host memory. */
36 if (location != OSKAR_CPU)
37 {
38 *status = OSKAR_ERR_BAD_LOCATION;
39 return;
40 }
41
42 if (type == OSKAR_SINGLE)
43 {
44 float *ra_ = 0, *dec_ = 0, *I_ = 0, *Q_ = 0, *U_ = 0, *V_ = 0;
45 float *ref_ = 0, *spix_ = 0, *rm_ = 0;
46 float *l_ = 0, *m_ = 0, *n_ = 0, *maj_ = 0, *min_ = 0, *pa_ = 0;
47 float *a_ = 0, *b_ = 0, *c_ = 0;
48 ra_ = oskar_mem_float(oskar_sky_ra_rad(sky), status);
49 dec_ = oskar_mem_float(oskar_sky_dec_rad(sky), status);
50 I_ = oskar_mem_float(oskar_sky_I(sky), status);
51 Q_ = oskar_mem_float(oskar_sky_Q(sky), status);
52 U_ = oskar_mem_float(oskar_sky_U(sky), status);
53 V_ = oskar_mem_float(oskar_sky_V(sky), status);
54 ref_ = oskar_mem_float(oskar_sky_reference_freq_hz(sky), status);
55 spix_ = oskar_mem_float(oskar_sky_spectral_index(sky), status);
56 rm_ = oskar_mem_float(oskar_sky_rotation_measure_rad(sky), status);
57 l_ = oskar_mem_float(oskar_sky_l(sky), status);
58 m_ = oskar_mem_float(oskar_sky_m(sky), status);
59 n_ = oskar_mem_float(oskar_sky_n(sky), status);
60 maj_ = oskar_mem_float(oskar_sky_fwhm_major_rad(sky), status);
61 min_ = oskar_mem_float(oskar_sky_fwhm_minor_rad(sky), status);
62 pa_ = oskar_mem_float(oskar_sky_position_angle_rad(sky), status);
63 a_ = oskar_mem_float(oskar_sky_gaussian_a(sky), status);
64 b_ = oskar_mem_float(oskar_sky_gaussian_b(sky), status);
65 c_ = oskar_mem_float(oskar_sky_gaussian_c(sky), status);
66
67 for (in = 0; in < num_sources; ++in)
68 {
69 if (!(I_[in] > (float)min_I && I_[in] <= (float)max_I)) continue;
70
71 ra_[out] = ra_[in];
72 dec_[out] = dec_[in];
73 I_[out] = I_[in];
74 Q_[out] = Q_[in];
75 U_[out] = U_[in];
76 V_[out] = V_[in];
77 ref_[out] = ref_[in];
78 spix_[out] = spix_[in];
79 rm_[out] = rm_[in];
80 l_[out] = l_[in];
81 m_[out] = m_[in];
82 n_[out] = n_[in];
83 maj_[out] = maj_[in];
84 min_[out] = min_[in];
85 pa_[out] = pa_[in];
86 a_[out] = a_[in];
87 b_[out] = b_[in];
88 c_[out] = c_[in];
89 out++;
90 }
91 }
92 else if (type == OSKAR_DOUBLE)
93 {
94 double *ra_ = 0, *dec_ = 0, *I_ = 0, *Q_ = 0, *U_ = 0, *V_ = 0;
95 double *ref_ = 0, *spix_ = 0, *rm_ = 0;
96 double *l_ = 0, *m_ = 0, *n_ = 0, *maj_ = 0, *min_ = 0, *pa_ = 0;
97 double *a_ = 0, *b_ = 0, *c_ = 0;
98 ra_ = oskar_mem_double(oskar_sky_ra_rad(sky), status);
99 dec_ = oskar_mem_double(oskar_sky_dec_rad(sky), status);
100 I_ = oskar_mem_double(oskar_sky_I(sky), status);
101 Q_ = oskar_mem_double(oskar_sky_Q(sky), status);
102 U_ = oskar_mem_double(oskar_sky_U(sky), status);
103 V_ = oskar_mem_double(oskar_sky_V(sky), status);
104 ref_ = oskar_mem_double(oskar_sky_reference_freq_hz(sky), status);
105 spix_ = oskar_mem_double(oskar_sky_spectral_index(sky), status);
106 rm_ = oskar_mem_double(oskar_sky_rotation_measure_rad(sky), status);
107 l_ = oskar_mem_double(oskar_sky_l(sky), status);
108 m_ = oskar_mem_double(oskar_sky_m(sky), status);
109 n_ = oskar_mem_double(oskar_sky_n(sky), status);
110 maj_ = oskar_mem_double(oskar_sky_fwhm_major_rad(sky), status);
111 min_ = oskar_mem_double(oskar_sky_fwhm_minor_rad(sky), status);
112 pa_ = oskar_mem_double(oskar_sky_position_angle_rad(sky), status);
113 a_ = oskar_mem_double(oskar_sky_gaussian_a(sky), status);
114 b_ = oskar_mem_double(oskar_sky_gaussian_b(sky), status);
115 c_ = oskar_mem_double(oskar_sky_gaussian_c(sky), status);
116
117 for (in = 0; in < num_sources; ++in)
118 {
119 if (!(I_[in] > min_I && I_[in] <= max_I)) continue;
120
121 ra_[out] = ra_[in];
122 dec_[out] = dec_[in];
123 I_[out] = I_[in];
124 Q_[out] = Q_[in];
125 U_[out] = U_[in];
126 V_[out] = V_[in];
127 ref_[out] = ref_[in];
128 spix_[out] = spix_[in];
129 rm_[out] = rm_[in];
130 l_[out] = l_[in];
131 m_[out] = m_[in];
132 n_[out] = n_[in];
133 maj_[out] = maj_[in];
134 min_[out] = min_[in];
135 pa_[out] = pa_[in];
136 a_[out] = a_[in];
137 b_[out] = b_[in];
138 c_[out] = c_[in];
139 out++;
140 }
141 }
142 else
143 {
144 *status = OSKAR_ERR_BAD_DATA_TYPE;
145 return;
146 }
147
148 /* Set the new size of the sky model. */
149 oskar_sky_resize(sky, out, status);
150 }
151
152 #ifdef __cplusplus
153 }
154 #endif
155