1 /*
2 * Copyright (c) 2016-2021, The OSKAR Developers.
3 * See the LICENSE file at the top-level directory of this distribution.
4 */
5
6 #include "imager/oskar_imager_linear_to_stokes.h"
7
8 #include <stddef.h>
9
10 #ifdef __cplusplus
11 extern "C" {
12 #endif
13
oskar_imager_linear_to_stokes(const oskar_Mem * in,oskar_Mem ** out,int * status)14 void oskar_imager_linear_to_stokes(const oskar_Mem* in, oskar_Mem** out,
15 int* status)
16 {
17 size_t i = 0;
18 if (*status) return;
19 const size_t num = oskar_mem_length(in);
20
21 /* Create output array or resize if required. */
22 if (!*out)
23 {
24 *out = oskar_mem_create(oskar_mem_type(in), OSKAR_CPU, num, status);
25 }
26 else
27 {
28 oskar_mem_ensure(*out, num, status);
29 }
30
31 /* Copy or convert if required. */
32 if (!oskar_mem_is_matrix(in))
33 {
34 /* Already Stokes I. */
35 oskar_mem_copy_contents(*out, in, 0, 0, num, status);
36 }
37 else
38 {
39 if (oskar_mem_precision(in) == OSKAR_DOUBLE)
40 {
41 const double4c* d_ = oskar_mem_double4c_const(in, status);
42 double4c* s_ = oskar_mem_double4c(*out, status);
43 for (i = 0; i < num; ++i)
44 {
45 /* I = 0.5 (XX + YY) */
46 s_[i].a.x = 0.5 * (d_[i].a.x + d_[i].d.x);
47 s_[i].a.y = 0.5 * (d_[i].a.y + d_[i].d.y);
48 /* Q = 0.5 (XX - YY) */
49 s_[i].b.x = 0.5 * (d_[i].a.x - d_[i].d.x);
50 s_[i].b.y = 0.5 * (d_[i].a.y - d_[i].d.y);
51 /* U = 0.5 (XY + YX) */
52 s_[i].c.x = 0.5 * (d_[i].b.x + d_[i].c.x);
53 s_[i].c.y = 0.5 * (d_[i].b.y + d_[i].c.y);
54 /* V = -0.5i (XY - YX) */
55 s_[i].d.x = 0.5 * (d_[i].b.y - d_[i].c.y);
56 s_[i].d.y = -0.5 * (d_[i].b.x - d_[i].c.x);
57 }
58 }
59 else
60 {
61 const float4c* d_ = oskar_mem_float4c_const(in, status);
62 float4c* s_ = oskar_mem_float4c(*out, status);
63 for (i = 0; i < num; ++i)
64 {
65 /* I = 0.5 (XX + YY) */
66 s_[i].a.x = 0.5 * (d_[i].a.x + d_[i].d.x);
67 s_[i].a.y = 0.5 * (d_[i].a.y + d_[i].d.y);
68 /* Q = 0.5 (XX - YY) */
69 s_[i].b.x = 0.5 * (d_[i].a.x - d_[i].d.x);
70 s_[i].b.y = 0.5 * (d_[i].a.y - d_[i].d.y);
71 /* U = 0.5 (XY + YX) */
72 s_[i].c.x = 0.5 * (d_[i].b.x + d_[i].c.x);
73 s_[i].c.y = 0.5 * (d_[i].b.y + d_[i].c.y);
74 /* V = -0.5i (XY - YX) */
75 s_[i].d.x = 0.5 * (d_[i].b.y - d_[i].c.y);
76 s_[i].d.y = -0.5 * (d_[i].b.x - d_[i].c.x);
77 }
78 }
79 }
80 }
81
82 #ifdef __cplusplus
83 }
84 #endif
85