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
8 #include "math/oskar_fit_ellipse.h"
9 #include "math/oskar_rotate.h"
10 #include "convert/oskar_convert_lon_lat_to_relative_directions.h"
11 #include "convert/oskar_convert_lon_lat_to_xyz.h"
12 #include "convert/oskar_convert_relative_directions_to_lon_lat.h"
13 #include "convert/oskar_convert_xyz_to_lon_lat.h"
14
15 #include <stdlib.h>
16 #include "math/oskar_cmath.h"
17
18 #define M_PI_2_2_LN_2 7.11941466249375271693034 /* pi^2 / (2 log_e(2)) */
19 #define MAX(x, y) (((x) > (y)) ? (x) : (y))
20 #define MIN(x, y) (((x) < (y)) ? (x) : (y))
21
22 #ifdef __cplusplus
23 extern "C" {
24 #endif
25
26 /* Number of points that define the ellipse */
27 #define ELLIPSE_PTS 6
28
oskar_sky_evaluate_gaussian_source_parameters(oskar_Sky * sky,int zero_failed_sources,double ra0,double dec0,int * num_failed,int * status)29 void oskar_sky_evaluate_gaussian_source_parameters(oskar_Sky* sky,
30 int zero_failed_sources, double ra0, double dec0, int* num_failed,
31 int* status)
32 {
33 int i = 0, j = 0;
34 if (*status) return;
35 if (oskar_sky_mem_location(sky) != OSKAR_CPU)
36 {
37 *status = OSKAR_ERR_BAD_LOCATION;
38 return;
39 }
40 const int type = oskar_sky_precision(sky);
41 const int num_sources = oskar_sky_num_sources(sky);
42 const double sin_dec0 = sin(dec0), cos_dec0 = cos(dec0);
43 if (type == OSKAR_DOUBLE)
44 {
45 /* Double precision. */
46 const double *ra_ = 0, *dec_ = 0, *maj_ = 0, *min_ = 0, *pa_ = 0;
47 double *I_ = 0, *Q_ = 0, *U_ = 0, *V_ = 0, *a_ = 0, *b_ = 0, *c_ = 0;
48 double cos_pa_2 = 0.0, sin_pa_2 = 0.0, sin_2pa = 0.0;
49 double inv_std_min_2 = 0.0, inv_std_maj_2 = 0.0;
50 double ellipse_a = 0.0, ellipse_b = 0.0, maj = 0.0, min = 0.0, pa = 0.0;
51 double cos_pa = 0.0, sin_pa = 0.0, t = 0.0;
52 double l[ELLIPSE_PTS], m[ELLIPSE_PTS];
53 double work1[5 * ELLIPSE_PTS], work2[5 * ELLIPSE_PTS];
54 double lon[ELLIPSE_PTS], lat[ELLIPSE_PTS];
55 double x[ELLIPSE_PTS], y[ELLIPSE_PTS], z[ELLIPSE_PTS];
56 ra_ = oskar_mem_double_const(oskar_sky_ra_rad_const(sky), status);
57 dec_ = oskar_mem_double_const(oskar_sky_dec_rad_const(sky), status);
58 maj_ = oskar_mem_double_const(oskar_sky_fwhm_major_rad_const(sky), status);
59 min_ = oskar_mem_double_const(oskar_sky_fwhm_minor_rad_const(sky), status);
60 pa_ = oskar_mem_double_const(oskar_sky_position_angle_rad_const(sky), status);
61 I_ = oskar_mem_double(oskar_sky_I(sky), status);
62 Q_ = oskar_mem_double(oskar_sky_Q(sky), status);
63 U_ = oskar_mem_double(oskar_sky_U(sky), status);
64 V_ = oskar_mem_double(oskar_sky_V(sky), status);
65 a_ = oskar_mem_double(oskar_sky_gaussian_a(sky), status);
66 b_ = oskar_mem_double(oskar_sky_gaussian_b(sky), status);
67 c_ = oskar_mem_double(oskar_sky_gaussian_c(sky), status);
68
69 for (i = 0; i < num_sources; ++i)
70 {
71 /* Note: could do something different from the projection below
72 * in the case of a line (i.e. maj or min = 0), as in this case
73 * there is no ellipse to project, only two points.
74 * -- This continue could then be a if() .. else() instead.
75 */
76 if (maj_[i] == 0.0 && min_[i] == 0.0) continue;
77
78 /* Evaluate shape of ellipse on the l,m plane. */
79 ellipse_a = maj_[i]/2.0;
80 ellipse_b = min_[i]/2.0;
81 cos_pa = cos(pa_[i]);
82 sin_pa = sin(pa_[i]);
83 #if 0
84 /* TODO(FD) Consider replacing existing code with this version? */
85 /* Get source parameters. */
86 const double cos_ra = cos(ra_[i]);
87 const double sin_ra = sin(ra_[i]);
88 const double cos_dec = cos(dec_[i]);
89 const double sin_dec = sin(dec_[i]);
90 for (j = 0; j < ELLIPSE_PTS; ++j)
91 {
92 double sin_t, cos_t;
93 /* Evaluate shape of ellipse on the l,m plane
94 * at RA = 0, Dec = 0. */
95 t = j * 60.0 * M_PI / 180.0; sin_t = sin(t); cos_t = cos(t);
96 const double l_ = ellipse_a * cos_t * sin_pa + ellipse_b * sin_t * cos_pa;
97 const double m_ = ellipse_a * cos_t * cos_pa - ellipse_b * sin_t * sin_pa;
98 const double n_ = sqrt(1.0 - l_*l_ - m_*m_);
99
100 /* Rotate to source position.
101 * For a source on tangent plane at (x, y, z) = (1, 0, 0),
102 * x parallel to n, y parallel to l, z parallel to m. */
103 const double x = n_ * cos_ra * cos_dec - l_ * sin_ra - m_ * cos_ra * sin_dec;
104 const double y = n_ * cos_dec * sin_ra + l_ * cos_ra - m_ * sin_ra * sin_dec;
105 const double z = n_ * sin_dec + m_ * cos_dec;
106
107 /* Get ellipse points relative to phase centre.
108 * sin(lat) is z; cos(lat) is sqrt(x*x + y*y); t is rel_lon */
109 t = atan2(y, x) - ra0; sin_t = sin(t); cos_t = cos(t);
110 const double cos_lat = sqrt(x*x + y*y);
111 l[j] = cos_lat * sin_t;
112 m[j] = cos_dec0 * z - sin_dec0 * cos_lat * cos_t;
113 }
114 #else
115 for (j = 0; j < ELLIPSE_PTS; ++j)
116 {
117 t = j * 60.0 * M_PI / 180.0;
118 l[j] = ellipse_a*cos(t)*sin_pa + ellipse_b*sin(t)*cos_pa;
119 m[j] = ellipse_a*cos(t)*cos_pa - ellipse_b*sin(t)*sin_pa;
120 }
121 oskar_convert_relative_directions_to_lon_lat_2d_d(ELLIPSE_PTS,
122 l, m, 0, 0.0, 1.0, 0.0, lon, lat);
123
124 /* Rotate on the sphere. */
125 oskar_convert_lon_lat_to_xyz_d(ELLIPSE_PTS, lon, lat, x, y, z);
126 oskar_rotate_sph_d(ELLIPSE_PTS, x, y, z, ra_[i], dec_[i]);
127 oskar_convert_xyz_to_lon_lat_d(ELLIPSE_PTS, x, y, z, lon, lat);
128
129 oskar_convert_lon_lat_to_relative_directions_2d_d(
130 ELLIPSE_PTS, lon, lat, ra0, cos_dec0, sin_dec0, l, m, 0);
131 #endif
132 /* Get new major and minor axes and position angle. */
133 oskar_fit_ellipse_d(&maj, &min, &pa, ELLIPSE_PTS, l, m, work1,
134 work2, status);
135
136 /* Check if fitting failed. */
137 if (*status == OSKAR_ERR_ELLIPSE_FIT_FAILED)
138 {
139 if (zero_failed_sources)
140 {
141 I_[i] = 0.0;
142 Q_[i] = 0.0;
143 U_[i] = 0.0;
144 V_[i] = 0.0;
145 }
146 ++(*num_failed);
147 *status = 0;
148 continue;
149 }
150 else if (*status) break;
151
152 /* Evaluate ellipse parameters. */
153 inv_std_maj_2 = 0.5 * (maj * maj) * M_PI_2_2_LN_2;
154 inv_std_min_2 = 0.5 * (min * min) * M_PI_2_2_LN_2;
155 cos_pa_2 = cos(pa) * cos(pa);
156 sin_pa_2 = sin(pa) * sin(pa);
157 sin_2pa = sin(2.0 * pa);
158 a_[i] = cos_pa_2*inv_std_min_2 + sin_pa_2*inv_std_maj_2;
159 b_[i] = -sin_2pa*inv_std_min_2*0.5 + sin_2pa *inv_std_maj_2*0.5;
160 c_[i] = sin_pa_2*inv_std_min_2 + cos_pa_2*inv_std_maj_2;
161 }
162 }
163 else
164 {
165 /* Single precision. */
166 const float *ra_ = 0, *dec_ = 0, *maj_ = 0, *min_ = 0, *pa_ = 0;
167 float *I_ = 0, *Q_ = 0, *U_ = 0, *V_ = 0, *a_ = 0, *b_ = 0, *c_ = 0;
168 float cos_pa_2 = 0.0, sin_pa_2 = 0.0, sin_2pa = 0.0;
169 float inv_std_min_2 = 0.0, inv_std_maj_2 = 0.0;
170 float ellipse_a = 0.0, ellipse_b = 0.0, maj = 0.0, min = 0.0, pa = 0.0;
171 float cos_pa = 0.0, sin_pa = 0.0, t = 0.0;
172 float l[ELLIPSE_PTS], m[ELLIPSE_PTS];
173 float work1[5 * ELLIPSE_PTS], work2[5 * ELLIPSE_PTS];
174 float lon[ELLIPSE_PTS], lat[ELLIPSE_PTS];
175 float x[ELLIPSE_PTS], y[ELLIPSE_PTS], z[ELLIPSE_PTS];
176 ra_ = oskar_mem_float_const(oskar_sky_ra_rad_const(sky), status);
177 dec_ = oskar_mem_float_const(oskar_sky_dec_rad_const(sky), status);
178 maj_ = oskar_mem_float_const(oskar_sky_fwhm_major_rad_const(sky), status);
179 min_ = oskar_mem_float_const(oskar_sky_fwhm_minor_rad_const(sky), status);
180 pa_ = oskar_mem_float_const(oskar_sky_position_angle_rad_const(sky), status);
181 I_ = oskar_mem_float(oskar_sky_I(sky), status);
182 Q_ = oskar_mem_float(oskar_sky_Q(sky), status);
183 U_ = oskar_mem_float(oskar_sky_U(sky), status);
184 V_ = oskar_mem_float(oskar_sky_V(sky), status);
185 a_ = oskar_mem_float(oskar_sky_gaussian_a(sky), status);
186 b_ = oskar_mem_float(oskar_sky_gaussian_b(sky), status);
187 c_ = oskar_mem_float(oskar_sky_gaussian_c(sky), status);
188
189 for (i = 0; i < num_sources; ++i)
190 {
191 /* Note: could do something different from the projection below
192 * in the case of a line (i.e. maj or min = 0), as in this case
193 * there is no ellipse to project, only two points.
194 * -- This continue could then be a if() .. else() instead.
195 */
196 if (maj_[i] == 0.0 && min_[i] == 0.0) continue;
197
198 /* Evaluate shape of ellipse on the l,m plane. */
199 ellipse_a = maj_[i]/2.0;
200 ellipse_b = min_[i]/2.0;
201 cos_pa = cos(pa_[i]);
202 sin_pa = sin(pa_[i]);
203 for (j = 0; j < ELLIPSE_PTS; ++j)
204 {
205 t = j * 60.0 * M_PI / 180.0;
206 l[j] = ellipse_a*cos(t)*sin_pa + ellipse_b*sin(t)*cos_pa;
207 m[j] = ellipse_a*cos(t)*cos_pa - ellipse_b*sin(t)*sin_pa;
208 }
209 oskar_convert_relative_directions_to_lon_lat_2d_f(ELLIPSE_PTS,
210 l, m, 0, 0.0, 1.0, 0.0, lon, lat);
211
212 /* Rotate on the sphere. */
213 oskar_convert_lon_lat_to_xyz_f(ELLIPSE_PTS, lon, lat, x, y, z);
214 oskar_rotate_sph_f(ELLIPSE_PTS, x, y, z, ra_[i], dec_[i]);
215 oskar_convert_xyz_to_lon_lat_f(ELLIPSE_PTS, x, y, z, lon, lat);
216
217 oskar_convert_lon_lat_to_relative_directions_2d_f(ELLIPSE_PTS,
218 lon, lat, (float)ra0, (float)cos_dec0, (float)sin_dec0,
219 l, m, 0);
220
221 /* Get new major and minor axes and position angle. */
222 oskar_fit_ellipse_f(&maj, &min, &pa, ELLIPSE_PTS, l, m, work1,
223 work2, status);
224
225 /* Check if fitting failed. */
226 if (*status == OSKAR_ERR_ELLIPSE_FIT_FAILED)
227 {
228 if (zero_failed_sources)
229 {
230 I_[i] = 0.0;
231 Q_[i] = 0.0;
232 U_[i] = 0.0;
233 V_[i] = 0.0;
234 }
235 ++(*num_failed);
236 *status = 0;
237 continue;
238 }
239 else if (*status) break;
240
241 /* Evaluate ellipse parameters. */
242 inv_std_maj_2 = 0.5 * (maj * maj) * M_PI_2_2_LN_2;
243 inv_std_min_2 = 0.5 * (min * min) * M_PI_2_2_LN_2;
244 cos_pa_2 = cos(pa) * cos(pa);
245 sin_pa_2 = sin(pa) * sin(pa);
246 sin_2pa = sin(2.0 * pa);
247 a_[i] = cos_pa_2*inv_std_min_2 + sin_pa_2*inv_std_maj_2;
248 b_[i] = -sin_2pa*inv_std_min_2*0.5 + sin_2pa *inv_std_maj_2*0.5;
249 c_[i] = sin_pa_2*inv_std_min_2 + cos_pa_2*inv_std_maj_2;
250 }
251 }
252 }
253
254 #ifdef __cplusplus
255 }
256 #endif
257