1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 #include <cdi.h>
9
10 #include "process_int.h"
11 #include "field.h"
12
13 static void
field2_add_complex(Field & field1,const Field & field2)14 field2_add_complex(Field &field1, const Field &field2)
15 {
16 const auto missval1 = field1.missval;
17 const auto missval2 = field2.missval;
18 auto &array1 = field1.vec_d;
19 const auto &array2 = field2.vec_d;
20
21 if (field1.nwpv != 2) cdo_abort("Field1 is not complex!");
22 if (field2.nwpv != 2) cdo_abort("Field2 is not complex!");
23
24 const auto gridsize = gridInqSize(field1.grid);
25 if (gridsize != gridInqSize(field2.grid)) cdo_abort("Fields have different size (%s)", __func__);
26
27 for (size_t i = 0; i < gridsize; i++)
28 {
29 array1[2 * i] = ADDMN(array1[2 * i], array2[2 * i]);
30 array1[2 * i + 1] = ADDMN(array1[2 * i + 1], array2[2 * i + 1]);
31 }
32 }
33
34 static void
field2_sub_complex(Field & field1,const Field & field2)35 field2_sub_complex(Field &field1, const Field &field2)
36 {
37 const auto missval1 = field1.missval;
38 const auto missval2 = field2.missval;
39 auto &array1 = field1.vec_d;
40 const auto &array2 = field2.vec_d;
41
42 if (field1.nwpv != 2) cdo_abort("Field1 is not complex!");
43 if (field2.nwpv != 2) cdo_abort("Field2 is not complex!");
44
45 const auto gridsize = gridInqSize(field1.grid);
46 if (gridsize != gridInqSize(field2.grid)) cdo_abort("Fields have different size (%s)", __func__);
47
48 for (size_t i = 0; i < gridsize; i++)
49 {
50 array1[2 * i] = SUBMN(array1[2 * i], array2[2 * i]);
51 array1[2 * i + 1] = SUBMN(array1[2 * i + 1], array2[2 * i + 1]);
52 }
53 }
54
55 static void
field2_mul_complex(Field & field1,const Field & field2)56 field2_mul_complex(Field &field1, const Field &field2)
57 {
58 const auto missval1 = field1.missval;
59 const auto missval2 = field2.missval;
60 auto &array1 = field1.vec_d;
61 const auto &array2 = field2.vec_d;
62
63 if (field1.nwpv != 2) cdo_abort("Field1 is not complex!");
64 if (field2.nwpv != 2) cdo_abort("Field2 is not complex!");
65
66 const auto gridsize = gridInqSize(field1.grid);
67 if (gridsize != gridInqSize(field2.grid)) cdo_abort("Fields have different size (%s)", __func__);
68
69 // z1 x z2 = (x1x2 - y1y2) + i(x1y2 + x2y1)
70 for (size_t i = 0; i < gridsize; i++)
71 {
72 const auto a1r = array1[2 * i];
73 const auto a1i = array1[2 * i + 1];
74 array1[2 * i] = SUBMN(MULMN(a1r, array2[2 * i]), MULMN(a1i, array2[2 * i + 1]));
75 array1[2 * i + 1] = ADDMN(MULMN(a1r, array2[2 * i + 1]), MULMN(a1i, array2[2 * i]));
76 }
77 }
78
79 static void
field2_div_complex(Field & field1,const Field & field2)80 field2_div_complex(Field &field1, const Field &field2)
81 {
82 const auto missval1 = field1.missval;
83 const auto missval2 = field2.missval;
84 auto &array1 = field1.vec_d;
85 const auto &array2 = field2.vec_d;
86
87 if (field1.nwpv != 2) cdo_abort("Field1 is not complex!");
88 if (field2.nwpv != 2) cdo_abort("Field2 is not complex!");
89
90 const auto gridsize = gridInqSize(field1.grid);
91 if (gridsize != gridInqSize(field2.grid)) cdo_abort("Fields have different size (%s)", __func__);
92
93 // z1 / z2 = (x1x2 + y1y2) / (x2x2 + y2y2) + i (y1x2 - x1y2) / (x2x2 + y2y2)
94 for (size_t i = 0; i < gridsize; i++)
95 {
96 const auto a1r = array1[2 * i];
97 const auto a1i = array1[2 * i + 1];
98 const auto denominator = ADDMN(MULMN(array2[2 * i], array2[2 * i]), MULMN(array2[2 * i + 1], array2[2 * i + 1]));
99 array1[2 * i] = DIVMN(ADDMN(MULMN(a1r, array2[2 * i]), MULMN(a1i, array2[2 * i + 1])), denominator);
100 array1[2 * i + 1] = DIVMN(SUBMN(MULMN(a1i, array2[2 * i]), MULMN(a1r, array2[2 * i + 1])), denominator);
101 }
102 }
103
104 void
field2_function_complex(Field & field1,const Field & field2,int function)105 field2_function_complex(Field &field1, const Field &field2, int function)
106 {
107 // clang-format off
108 switch (function)
109 {
110 case FieldFunc_Add: field2_add_complex(field1, field2); break;
111 case FieldFunc_Sub: field2_sub_complex(field1, field2); break;
112 case FieldFunc_Mul: field2_mul_complex(field1, field2); break;
113 case FieldFunc_Div: field2_div_complex(field1, field2); break;
114 default: cdo_abort("%s: function %d not implemented!", __func__, function);
115 }
116 // clang-format on
117 }
118