YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_lat_clipping.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <stdlib.h>
6#include <stdio.h>
7#include <math.h>
8
9#include "grids/basic_grid.h"
10#include "clipping.h"
11#include "geometry.h"
12#include "tests.h"
13#include "area.h"
14#include "test_common.h"
15
21static void utest_lat_clipping(struct yac_grid_cell cells,
22 double deg_lat_bounds[2],
23 struct yac_grid_cell * ref_cells,
24 unsigned num_ref_cells);
25
26static int compare_cells(struct yac_grid_cell a, struct yac_grid_cell b);
27
28static double utest_get_intersection(
29 enum yac_edge_type edge_type, double points[2][2], double lat) {
30
31 double intersection[3], intersection_lon, intersection_lat;
32
33 if (!intersect(edge_type, points[0][0], points[0][1],
34 points[1][0], points[1][1],
35 YAC_LAT_CIRCLE_EDGE, points[0][0], lat,
36 points[1][0], lat, intersection))
37 exit(EXIT_FAILURE);
38 XYZtoLL(intersection, &intersection_lon, &intersection_lat);
39 return intersection_lon / YAC_RAD;
40}
41
42int main (void) {
43
44 enum yac_edge_type gc_edges[] = {
45 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
46 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
47 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
48 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE};
49 enum yac_edge_type latlon_edges[] = {
50 YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
51
52 double temp_lon[16];
53
54 // a gaussian cell outside the lat bound
55 utest_lat_clipping(
56 generate_cell_deg(
57 (double[4]){-5,5,5,-5},
58 (double[4]){-5,-5,5,5}, latlon_edges, 4),
59 (double[2]){10, 20}, NULL, 0);
60
61 // a gaussian cell outside the lat bound
62 utest_lat_clipping(
63 generate_cell_deg(
64 (double[4]){-5,5,5,-5},
65 (double[4]){-5,-5,10,10}, latlon_edges, 4),
66 (double[2]){10, 20}, NULL, 0);
67
68 // a gaussian cell outside the lat bound
69 utest_lat_clipping(
70 generate_cell_deg(
71 (double[4]){-5,5,5,-5},
72 (double[4]){20,20,30,30}, latlon_edges, 4),
73 (double[2]){10, 20}, NULL, 0);
74
75 // a gaussian cell outside the lat bound
76 utest_lat_clipping(
77 generate_cell_deg(
78 (double[4]){-5,5,5,-5},
79 (double[4]){80,80,90,90}, latlon_edges, 4),
80 (double[2]){10, 20}, NULL, 0);
81
82 // a gaussian cell outside the lat bound
83 utest_lat_clipping(
84 generate_cell_deg(
85 (double[4]){-5,5,5,-5},
86 (double[4]){-80,-80,-90,-90}, latlon_edges, 4),
87 (double[2]){10, 20}, NULL, 0);
88
89 // a gaussian cell outside the lat bound (one bound is at a pole)
90 utest_lat_clipping(
91 generate_cell_deg(
92 (double[4]){-5,5,5,-5},
93 (double[4]){10,10,20,20}, latlon_edges, 4),
94 (double[2]){80, 90}, NULL, 0);
95
96 // a gaussian cell outside the lat bound (both bounds are identical)
97 utest_lat_clipping(
98 generate_cell_deg(
99 (double[4]){-5,5,5,-5},
100 (double[4]){10,10,20,20}, latlon_edges, 4),
101 (double[2]){0, 0}, NULL, 0);
102
103 // a gaussian cell outside the lat bound (both bounds are identical)
104 utest_lat_clipping(
105 generate_cell_deg(
106 (double[4]){-5,5,5,-5},
107 (double[4]){-10,-10,-20,-20}, latlon_edges, 4),
108 (double[2]){0, 0}, NULL, 0);
109
110 // a gaussian cell outside the lat bound (both bounds are identical)
111 utest_lat_clipping(
112 generate_cell_deg(
113 (double[4]){-5,5,5,-5},
114 (double[4]){-5,-5,5,5}, latlon_edges, 4),
115 (double[2]){0, 0}, NULL, 0);
116
117 // a gaussian cell overlapping with the lat bound
118 utest_lat_clipping(
119 generate_cell_deg(
120 (double[4]){-5,5,5,-5},
121 (double[4]){-10,-10,10,10}, latlon_edges, 4),
122 (double[2]){-5, 5},
123 (struct yac_grid_cell[]){
124 generate_cell_deg(
125 (double[4]){-5,5,5,-5},
126 (double[4]){-5,-5,5,5}, latlon_edges, 4)}, 1);
127
128 // a gaussian cell overlapping with the lat bound
129 utest_lat_clipping(
130 generate_cell_deg(
131 (double[4]){40,45,45,40},
132 (double[4]){10,10,20,20}, latlon_edges, 4),
133 (double[2]){12, 15},
134 (struct yac_grid_cell[]){
135 generate_cell_deg(
136 (double[4]){40,45,45,40},
137 (double[4]){12,12,15,15}, latlon_edges, 4)}, 1);
138
139 // a gaussian cell overlapping with the lat bound
140 utest_lat_clipping(
141 generate_cell_deg(
142 (double[4]){40,45,45,40},
143 (double[4]){10,10,20,20}, latlon_edges, 4),
144 (double[2]){15, 25},
145 (struct yac_grid_cell[]){
146 generate_cell_deg(
147 (double[4]){40,45,45,40},
148 (double[4]){15,15,20,20}, latlon_edges, 4)}, 1);
149
150 // a gaussian cell overlapping with the lat bound
151 utest_lat_clipping(
152 generate_cell_deg(
153 (double[4]){40,45,45,40},
154 (double[4]){10,10,20,20}, latlon_edges, 4),
155 (double[2]){15, 5},
156 (struct yac_grid_cell[]){
157 generate_cell_deg(
158 (double[4]){40,45,45,40},
159 (double[4]){15,15,10,10}, latlon_edges, 4)}, 1);
160
161 // a gaussian cell overlapping with the lat bound
162 utest_lat_clipping(
163 generate_cell_deg(
164 (double[4]){40,45,45,40},
165 (double[4]){10,10,20,20}, latlon_edges, 4),
166 (double[2]){5, 25},
167 (struct yac_grid_cell[]){
168 generate_cell_deg(
169 (double[4]){40,45,45,40},
170 (double[4]){10,10,20,20}, latlon_edges, 4)}, 1);
171
172 // a gaussian cell
173 utest_lat_clipping(
174 generate_cell_deg(
175 (double[4]){40,45,45,40},
176 (double[4]){10,10,20,20}, latlon_edges, 4),
177 (double[2]){90, 80}, NULL, 0);
178
179 // a gaussian cell
180 utest_lat_clipping(
181 generate_cell_deg(
182 (double[4]){40,45,45,40},
183 (double[4]){90,90,70,70}, latlon_edges, 4),
184 (double[2]){90, 80},
185 (struct yac_grid_cell[]){
186 generate_cell_deg(
187 (double[4]){40,45,45,40},
188 (double[4]){90,90,80,80}, latlon_edges, 4)}, 1);
189
190 // a gaussian cell
191 utest_lat_clipping(
192 generate_cell_deg(
193 (double[4]){40,45,45,40},
194 (double[4]){90,90,85,85}, latlon_edges, 4),
195 (double[2]){90, 80},
196 (struct yac_grid_cell[]){
197 generate_cell_deg(
198 (double[4]){40,45,45,40},
199 (double[4]){90,90,85,85}, latlon_edges, 4)}, 1);
200
201 // a gaussian cell
202 utest_lat_clipping(
203 generate_cell_deg(
204 (double[4]){40,45,45,40},
205 (double[4]){88,88,83,83}, latlon_edges, 4),
206 (double[2]){90, 80},
207 (struct yac_grid_cell[]){
208 generate_cell_deg(
209 (double[4]){40,45,45,40},
210 (double[4]){88,88,83,83}, latlon_edges, 4)}, 1);
211
212 // a gaussian cell
213 utest_lat_clipping(
214 generate_cell_deg(
215 (double[4]){40,45,45,40},
216 (double[4]){88,88,83,83}, latlon_edges, 4),
217 (double[2]){90, 80},
218 (struct yac_grid_cell[]){
219 generate_cell_deg(
220 (double[4]){40,45,45,40},
221 (double[4]){88,88,83,83}, latlon_edges, 4)}, 1);
222
223 // a gaussian cell
224 utest_lat_clipping(
225 generate_cell_deg(
226 (double[4]){40,45,45,40},
227 (double[4]){88,88,70,70}, latlon_edges, 4),
228 (double[2]){90, 80},
229 (struct yac_grid_cell[]){
230 generate_cell_deg(
231 (double[4]){40,45,45,40},
232 (double[4]){88,88,80,80}, latlon_edges, 4)}, 1);
233
234 // a gaussian cell
235 utest_lat_clipping(
236 generate_cell_deg(
237 (double[4]){40,45,45,40},
238 (double[4]){10,10,20,20}, latlon_edges, 4),
239 (double[2]){-90, -80}, NULL, 0);
240
241 // a gaussian cell
242 utest_lat_clipping(
243 generate_cell_deg(
244 (double[4]){40,45,45,40},
245 (double[4]){-90,-90,-70,-70}, latlon_edges, 4),
246 (double[2]){-90, -80},
247 (struct yac_grid_cell[]){
248 generate_cell_deg(
249 (double[4]){40,45,45,40},
250 (double[4]){-90,-90,-80,-80}, latlon_edges, 4)}, 1);
251
252 // a gaussian cell
253 utest_lat_clipping(
254 generate_cell_deg(
255 (double[4]){40,45,45,40},
256 (double[4]){-90,-90,-85,-85}, latlon_edges, 4),
257 (double[2]){-90, -80},
258 (struct yac_grid_cell[]){
259 generate_cell_deg(
260 (double[4]){40,45,45,40},
261 (double[4]){-90,-90,-85,-85}, latlon_edges, 4)}, 1);
262
263 // a gaussian cell
264 utest_lat_clipping(
265 generate_cell_deg(
266 (double[4]){40,45,45,40},
267 (double[4]){-88,-88,-83,-83}, latlon_edges, 4),
268 (double[2]){-90, -80},
269 (struct yac_grid_cell[]){
270 generate_cell_deg(
271 (double[4]){40,45,45,40},
272 (double[4]){-88,-88,-83,-83}, latlon_edges, 4)}, 1);
273
274 // a gaussian cell
275 utest_lat_clipping(
276 generate_cell_deg(
277 (double[4]){40,45,45,40},
278 (double[4]){-88,-88,-83,-83}, latlon_edges, 4),
279 (double[2]){-90, -80},
280 (struct yac_grid_cell[]){
281 generate_cell_deg(
282 (double[4]){40,45,45,40},
283 (double[4]){-88,-88,-83,-83}, latlon_edges, 4)}, 1);
284
285 // a gaussian cell
286 utest_lat_clipping(
287 generate_cell_deg(
288 (double[4]){40,45,45,40},
289 (double[4]){-88,-88,-70,-70}, latlon_edges, 4),
290 (double[2]){-90, -80},
291 (struct yac_grid_cell[]){
292 generate_cell_deg(
293 (double[4]){40,45,45,40},
294 (double[4]){-88,-88,-80,-80}, latlon_edges, 4)}, 1);
295
296 // a gaussian cell that touches the north pole and partially overlaps
297 // with the bounds
298 utest_lat_clipping(
299 generate_cell_deg(
300 (double[4]){0,5,5,0},
301 (double[4]){90,90,85,85}, latlon_edges, 4),
302 (double[2]){89.9, 80},
303 (struct yac_grid_cell[]){
304 generate_cell_deg(
305 (double[4]){0,5,5,0},
306 (double[4]){89.9,89.9,85,85}, latlon_edges, 4)}, 1);
307
308 // a triangle cell
309 utest_lat_clipping(
310 generate_cell_deg(
311 (double[3]){10,20,15},
312 (double[3]){30,30,40}, gc_edges, 3),
313 (double[2]){35, 50},
314 (struct yac_grid_cell[]){
315 generate_cell_deg(
316 (double[]){15,
317 utest_get_intersection(
318 YAC_GREAT_CIRCLE_EDGE,
319 (double[2][2]){{15, 40}, {10, 30}}, 35),
320 utest_get_intersection(
321 YAC_GREAT_CIRCLE_EDGE,
322 (double[2][2]){{15, 40}, {20, 30}}, 35)},
323 (double[]){40,35,35},
324 (enum yac_edge_type[]){
325 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 3)}, 1);
326
327 // a triangle cell
328 utest_lat_clipping(
329 generate_cell_deg(
330 (double[3]){10,20,15},
331 (double[3]){30,30,40}, gc_edges, 3),
332 (double[2]){35, 25},
333 (struct yac_grid_cell[]){
334 generate_cell_deg(
335 (double[]){20, 10,
336 utest_get_intersection(
337 YAC_GREAT_CIRCLE_EDGE,
338 (double[2][2]){{15, 40}, {10, 30}}, 35),
339 utest_get_intersection(
340 YAC_GREAT_CIRCLE_EDGE,
341 (double[2][2]){{15, 40}, {20, 30}}, 35)},
342 (double[]){30,30,35,35},
343 (enum yac_edge_type[]){YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
344 YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 4)}, 1);
345
346 // a triangle cell
347 utest_lat_clipping(
348 generate_cell_deg(
349 (double[]){10,20,15},
350 (double[]){30,30,40}, gc_edges, 3),
351 (double[2]){45, 25},
352 (struct yac_grid_cell[]){
353 generate_cell_deg(
354 (double[]){10,20,15},
355 (double[]){30,30,40}, gc_edges, 3)}, 1);
356
357 // a triangle cell
358 utest_lat_clipping(
359 generate_cell_deg(
360 (double[3]){10,20,15},
361 (double[3]){30,30,40}, gc_edges, 3),
362 (double[2]){32, 38},
363 (struct yac_grid_cell[]){
364 generate_cell_deg(
365 (double[]){utest_get_intersection(
366 YAC_GREAT_CIRCLE_EDGE,
367 (double[2][2]){{15, 40}, {20, 30}}, 32),
368 utest_get_intersection(
369 YAC_GREAT_CIRCLE_EDGE,
370 (double[2][2]){{15, 40}, {20, 30}}, 38),
371 utest_get_intersection(
372 YAC_GREAT_CIRCLE_EDGE,
373 (double[2][2]){{15, 40}, {10, 30}}, 38),
374 utest_get_intersection(
375 YAC_GREAT_CIRCLE_EDGE,
376 (double[2][2]){{15, 40}, {10, 30}}, 32)},
377 (double[]){32,38,38,32},
378 (enum yac_edge_type[]){YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE,
379 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE}, 4)}, 1);
380
381 // a triangle cell
382 utest_lat_clipping(
383 generate_cell_deg(
384 (double[3]){-20,20,0},
385 (double[3]){80,80,85}, gc_edges, 3),
386 (double[2]){80.001, 89},
387 (struct yac_grid_cell[]){
388 generate_cell_deg(
389 (double[]){0,
390 utest_get_intersection(
391 YAC_GREAT_CIRCLE_EDGE,
392 (double[2][2]){{0, 85}, {-20, 80}}, 80.001),
393 -fabs(utest_get_intersection(
394 YAC_GREAT_CIRCLE_EDGE,
395 (double[2][2]){{-20, 80}, {20, 80}}, 80.001)),
396 fabs(utest_get_intersection(
397 YAC_GREAT_CIRCLE_EDGE,
398 (double[2][2]){{-20, 80}, {20, 80}}, 80.001)),
399 utest_get_intersection(
400 YAC_GREAT_CIRCLE_EDGE,
401 (double[2][2]){{0, 85}, {20, 80}}, 80.001),},
402 (double[]){85,80.001,80.001,80.001,80.001},
403 (enum yac_edge_type[]){YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE,
404 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE,
405 YAC_GREAT_CIRCLE_EDGE}, 5)}, 1);
406
407 // a triangle cell
408 utest_lat_clipping(
409 generate_cell_deg(
410 (double[]){10,20,15},
411 (double[]){82,82,88}, gc_edges, 3),
412 (double[2]){90, 80},
413 (struct yac_grid_cell[]){
414 generate_cell_deg(
415 (double[]){10,20,15},
416 (double[]){82,82,88}, gc_edges, 3)}, 1);
417
418 // a triangle cell
419 utest_lat_clipping(
420 generate_cell_deg(
421 (double[]){10,20,15},
422 (double[]){-82,-82,-88}, gc_edges, 3),
423 (double[2]){90, 80}, NULL, 0);
424
425 // a triangle cell
426 utest_lat_clipping(
427 generate_cell_deg(
428 (double[]){10,20,15},
429 (double[]){-82,-82,-88}, gc_edges, 3),
430 (double[2]){-90, -80},
431 (struct yac_grid_cell[]){
432 generate_cell_deg(
433 (double[]){10,20,15},
434 (double[]){-82,-82,-88}, gc_edges, 3)}, 1);
435
436 // a triangle cell
437 utest_lat_clipping(
438 generate_cell_deg(
439 (double[]){10,20,15},
440 (double[]){82,82,88}, gc_edges, 3),
441 (double[2]){-90, -80}, NULL, 0);
442
443 // a triangle cell
444 utest_lat_clipping(
445 generate_cell_deg(
446 (double[]){10,20,15},
447 (double[]){82,82,88}, gc_edges, 3),
448 (double[2]){90, 90}, NULL, 0);
449
450 // a triangle cell
451 utest_lat_clipping(
452 generate_cell_deg(
453 (double[]){10,20,15},
454 (double[]){82,82,88}, gc_edges, 3),
455 (double[2]){-90, -90}, NULL, 0);
456
457 // a triangle cell covering the pole
458 utest_lat_clipping(
459 generate_cell_deg(
460 (double[]){0,120,240},
461 (double[]){80,80,80}, gc_edges, 3),
462 (double[2]){90, 70},
463 (struct yac_grid_cell[]){
464 generate_cell_deg(
465 (double[]){0,120,240},
466 (double[]){80,80,80}, gc_edges, 3)}, 1);
467
468 // a triangle cell covering the pole
469 utest_lat_clipping(
470 generate_cell_deg(
471 (double[]){0,120,240},
472 (double[]){80,80,80}, gc_edges, 3),
473 (double[2]){60, 70}, NULL, 0);
474
475 // a triangle cell covering the pole
476 temp_lon[0] =
477 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
478 temp_lon[1] =
479 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
480 utest_lat_clipping(
481 generate_cell_deg(
482 (double[]){0,120,240},
483 (double[]){85,70,70}, gc_edges, 3),
484 (double[2]){90, 80},
485 (struct yac_grid_cell[]){
486 generate_cell_deg(
487 (double[]){0,temp_lon[0],temp_lon[1]},
488 (double[]){85,80,80},
489 (enum yac_edge_type[]){
490 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 3)}, 1);
491
492 // a triangle cell covering the pole
493 temp_lon[0] =
494 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 72);
495 temp_lon[1] =
496 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{120, 70}, {240, 70}}, 72);
497 temp_lon[1] = (temp_lon[1] < 0)?(temp_lon[1] + 360):(temp_lon[1]);
498 temp_lon[1] = (temp_lon[1] < 180)?(temp_lon[1]):(360 - temp_lon[1]);
499 temp_lon[2] = 360 - temp_lon[1];
500 temp_lon[3] =
501 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 72);
502 utest_lat_clipping(
503 generate_cell_deg(
504 (double[]){0,120,240},
505 (double[]){85,70,70}, gc_edges, 3),
506 (double[2]){90, 72},
507 (struct yac_grid_cell[]){
508 generate_cell_deg(
509 (double[]){0,temp_lon[0],temp_lon[1],temp_lon[2],temp_lon[3]},
510 (double[]){85,72,72,72,72},
511 (enum yac_edge_type[]){
512 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
513 YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 5)}, 1);
514
515 // a triangle cell covering the pole
516 temp_lon[0] =
517 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
518 temp_lon[1] =
519 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
520 utest_lat_clipping(
521 generate_cell_deg(
522 (double[]){0,120,240},
523 (double[]){85,70,70}, gc_edges, 3),
524 (double[2]){80, 65},
525 (struct yac_grid_cell[]){
526 generate_cell_deg(
527 (double[]){temp_lon[0],120,240,temp_lon[1]},
528 (double[]){80,70,70,80},
529 (enum yac_edge_type[]){
530 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
531 YAC_LAT_CIRCLE_EDGE}, 4)}, 1);
532
533 // a triangle cell covering the pole
534 temp_lon[0] =
535 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 84);
536 temp_lon[1] =
537 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
538 temp_lon[2] =
539 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
540 temp_lon[3] =
541 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 84);
542 utest_lat_clipping(
543 generate_cell_deg(
544 (double[]){0,120,240},
545 (double[]){85,70,70}, gc_edges, 3),
546 (double[2]){80, 84},
547 (struct yac_grid_cell[]){
548 generate_cell_deg(
549 (double[]){temp_lon[0],temp_lon[1],temp_lon[2],temp_lon[3]},
550 (double[]){84,80,80,84},
551 (enum yac_edge_type[]){
552 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
553 YAC_LAT_CIRCLE_EDGE}, 4)}, 1);
554
555 // a triangle cell covering the pole
556 temp_lon[0] =
557 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
558 temp_lon[1] =
559 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 72);
560 temp_lon[2] =
561 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{240, 70}, {120, 70}}, 72);
562 temp_lon[2] = (temp_lon[2] < 0)?(temp_lon[2] + 360):(temp_lon[2]);
563 temp_lon[2] = (temp_lon[2] < 180)?(temp_lon[2]):(360 - temp_lon[2]);
564 temp_lon[3] = 360 - temp_lon[2];
565 temp_lon[4] =
566 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 72);
567 temp_lon[5] =
568 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
569 utest_lat_clipping(
570 generate_cell_deg(
571 (double[]){0,120,240},
572 (double[]){85,70,70}, gc_edges, 3),
573 (double[2]){80, 72},
574 (struct yac_grid_cell[]){
575 generate_cell_deg(
576 (double[]){temp_lon[0],temp_lon[1],temp_lon[2],
577 temp_lon[3],temp_lon[4],temp_lon[5]},
578 (double[]){80,72,72,72,72,80},
579 (enum yac_edge_type[]){
580 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
581 YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE}, 6)}, 1);
582
583 // a triangle cell touching the pole
584 temp_lon[0] =
585 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 90}, {120, 80}}, 85);
586 temp_lon[1] =
587 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 90}, {240, 80}}, 85);
588 utest_lat_clipping(
589 generate_cell_deg(
590 (double[]){0,120,240},
591 (double[]){90,80,80}, gc_edges, 3),
592 (double[2]){90, 85},
593 (struct yac_grid_cell[]){
594 generate_cell_deg(
595 (double[]){0,temp_lon[0],temp_lon[1]},
596 (double[]){90,85,85},
597 (enum yac_edge_type[]){
598 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 3)}, 1);
599
600 // a triangle cell touching the pole
601 temp_lon[0] =
602 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 0}, {10, 5}}, 2.5);
603 temp_lon[1] =
604 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{10, -5}, {10, 5}}, 2.5);
605 temp_lon[2] =
606 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{10, -5}, {10, 5}}, -2.5);
607 temp_lon[3] =
608 utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{10, -5}, {0, 0}}, -2.5);
609 utest_lat_clipping(
610 generate_cell_deg(
611 (double[]){0,10,10},
612 (double[]){0,5,-5}, gc_edges, 3),
613 (double[2]){2.5, -2.5},
614 (struct yac_grid_cell[]){
615 generate_cell_deg(
616 (double[]){0,temp_lon[0],temp_lon[1],temp_lon[2],temp_lon[3]},
617 (double[]){0,2.5,2.5,-2.5,-2.5},
618 (enum yac_edge_type[]){
619 YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
620 YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 5)}, 1);
621
622 // empty great circle cell
623 utest_lat_clipping(
624 generate_cell_deg(
625 (double[]){0,0,0},
626 (double[]){-5,0,5}, gc_edges, 3),
627 (double[2]){2.5, -2.5}, NULL, 0);
628
629 return TEST_EXIT_CODE;
630}
631
632static int compare_cells(struct yac_grid_cell a, struct yac_grid_cell b) {
633
634 // Trigonometric functions can be very inaccurate depending on the compiler
635 // and the optimsation level. Therefore, we use angle_tol instead of
636 // yac_angle_tol in this routine.
637 double const angle_tol = 1e-6;
638
639 if ((yac_huiliers_area(a) < angle_tol * angle_tol) &&
640 (yac_huiliers_area(b) < angle_tol * angle_tol)) return 0;
641
642 if (a.num_corners != b.num_corners) return 1;
643
644 if (a.num_corners == 0) return 0;
645
646 for (int order = -1; order <= 1; order += 2) {
647 for (int start = 0; start < (int)(a.num_corners); ++start) {
648
649 int differences = 0;
650
651 for (int i = 0; i < (int)(a.num_corners); ++i) {
652
653 int j =
654 ((int)(a.num_corners) + start + order * i) % (int)(a.num_corners);
655
656 if (get_vector_angle(
657 a.coordinates_xyz[i], b.coordinates_xyz[j]) > angle_tol)
658 ++differences;
659
660 else if ((a.edge_type[i] == YAC_LAT_CIRCLE_EDGE) !=
661 (b.edge_type[
662 (((int)(a.num_corners))+j-(order<0))%
663 ((int)(a.num_corners))] == YAC_LAT_CIRCLE_EDGE))
664 ++differences;
665 }
666
667 if (!differences) return 0;
668 }
669 }
670
671 return 1;
672}
673
674static void utest_lat_clipping(struct yac_grid_cell cell,
675 double deg_lat_bounds[2],
676 struct yac_grid_cell * ref_cells,
677 unsigned num_ref_cells) {
678
679 double const area_tol = 1e-8 * yac_huiliers_area(cell);
680
681 double rad_lat_bounds[2] =
682 {deg_lat_bounds[0] * YAC_RAD, deg_lat_bounds[1] * YAC_RAD};
683
684 double mem_dummy_[128][3];
685 enum yac_edge_type edge_dummy[128];
686 struct yac_grid_cell overlap_cell,
687 test_cell =
688 (struct yac_grid_cell){.coordinates_xyz = mem_dummy_,
689 .edge_type = edge_dummy,
690 .num_corners = cell.num_corners};
691
692 yac_init_grid_cell(&overlap_cell);
693
694 for (int pole = 0; pole < 2; ++pole) {
695 for (int order = -1; order <= 1; order += 2) {
696 for (int start = 0; start < (int)(cell.num_corners); ++start) {
697 for (int i = 0; i < ((int)(cell.num_corners)); ++i) {
698 int j = (((int)(cell.num_corners))+start+i*order)%
699 ((int)(cell.num_corners));
700 test_cell.coordinates_xyz[i][0] = cell.coordinates_xyz[j][0];
701 test_cell.coordinates_xyz[i][1] = cell.coordinates_xyz[j][1];
702 test_cell.coordinates_xyz[i][2] = cell.coordinates_xyz[j][2];
703 j = (((int)(cell.num_corners)) + j - (order < 0))%
704 ((int)(cell.num_corners));
705 test_cell.edge_type[i] = cell.edge_type[j];
706 }
707
708 for (int bound_tol_a = -1; bound_tol_a <= 1; ++bound_tol_a) {
709 for (int bound_tol_b = -1; bound_tol_b <= 1; ++bound_tol_b) {
710 for (int bound_order = 0; bound_order <= 1; ++bound_order) {
711 double test_rad_lat_bounds[2] = {
712 MIN(M_PI_2,
713 MAX(-M_PI_2,
714 rad_lat_bounds[bound_order] +
715 (double)bound_tol_a * yac_angle_low_tol)),
716 MIN(M_PI_2,
717 MAX(-M_PI_2,
718 rad_lat_bounds[bound_order^1] +
719 (double)bound_tol_b * yac_angle_low_tol))};
720 yac_cell_lat_clipping(
721 1, &test_cell,test_rad_lat_bounds, &overlap_cell);
722 if (num_ref_cells == 0) {
723 if ((overlap_cell.num_corners != 0) &&
724 (yac_huiliers_area(overlap_cell) > area_tol))
725 PUT_ERR("ERROR: wrong clipping cell\n");
726 } else {
727 int match = 0;
728 for (int i = 0; i < (int)num_ref_cells; ++i)
729 match |= !compare_cells(overlap_cell, ref_cells[i]);
730 if (!match)
731 PUT_ERR("ERROR: wrong clipping cell\n");
732 }
733 }
734 }
735 }
736 }
737 }
738
739 // switch to the other pole
740 for (int i = 0; i < (int)(cell.num_corners); ++i)
741 cell.coordinates_xyz[i][2] *= -1.0;
742 rad_lat_bounds[0] *= -1.0, rad_lat_bounds[1] *= -1.0;
743 for (unsigned i = 0; i < num_ref_cells; ++i)
744 for (int j = 0; j < (int)(ref_cells[i].num_corners); ++j)
745 ref_cells[i].coordinates_xyz[j][2] *= -1.0;
746 }
747
748 yac_free_grid_cell(&cell);
749 for (unsigned i = 0; i < num_ref_cells; ++i)
750 yac_free_grid_cell(&ref_cells[i]);
751 yac_free_grid_cell(&overlap_cell);
752}
yac_huiliers_area
double yac_huiliers_area(struct yac_grid_cell cell)
Area calculation on a unit sphere taken from ESMF based on L'Huilier's Theorem.
Definition area.c:449
area.h
Structs and interfaces for area calculations.
basic_grid.h
yac_cell_lat_clipping
void yac_cell_lat_clipping(size_t N, struct yac_grid_cell *cells, double lat_bounds[2], struct yac_grid_cell *overlap_buffer)
cell clipping to get the cells describing the intersections
Definition clipping.c:1325
YAC_RAD
#define YAC_RAD
Definition dynamic_config.c:12
yac_angle_low_tol
#define yac_angle_low_tol
Definition geometry.h:29
get_vector_angle
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:340
yac_init_grid_cell
void yac_init_grid_cell(struct yac_grid_cell *cell)
Definition grid_cell.c:14
yac_free_grid_cell
void yac_free_grid_cell(struct yac_grid_cell *cell)
Definition grid_cell.c:44
yac_edge_type
yac_edge_type
Definition grid_cell.h:12
YAC_GREAT_CIRCLE_EDGE
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
YAC_LAT_CIRCLE_EDGE
@ YAC_LAT_CIRCLE_EDGE
latitude circle
Definition grid_cell.h:14
YAC_LON_CIRCLE_EDGE
@ YAC_LON_CIRCLE_EDGE
longitude circle
Definition grid_cell.h:15
test_raw_exchange.i
int i
Definition test_raw_exchange.py:89
XYZtoLL
static void XYZtoLL(double const p_in[], double *lon, double *lat)
Definition perf_toy_icon.c:355
yac_grid_cell
Definition grid_cell.h:18
yac_grid_cell::num_corners
size_t num_corners
Definition grid_cell.h:21
yac_grid_cell::edge_type
enum yac_edge_type * edge_type
Definition grid_cell.h:20
yac_grid_cell::coordinates_xyz
double(* coordinates_xyz)[3]
Definition grid_cell.h:19
generate_cell_deg
struct yac_grid_cell generate_cell_deg(double *lon, double *lat, enum yac_edge_type *edge_type, size_t num_corners)
Definition test_common.c:36
intersect
int intersect(enum yac_edge_type edge_type_a, double lon_a, double lat_a, double lon_b, double lat_b, enum yac_edge_type edge_type_b, double lon_c, double lat_c, double lon_d, double lat_d, double *intersection)
Definition test_common.c:64
test_common.h
gc_edges
static enum yac_edge_type gc_edges[]
Definition test_compute_overlap_area.c:16
latlon_edges
static enum yac_edge_type latlon_edges[4]
Definition test_compute_overlap_area.c:26
compare_cells
static int compare_cells(struct yac_grid_cell a, struct yac_grid_cell b)
Definition test_lat_clipping.c:632
TEST_EXIT_CODE
#define TEST_EXIT_CODE
Definition tests.h:14
PUT_ERR
#define PUT_ERR(string)
Definition tests.h:10
MIN
#define MIN(a, b)
Definition toy_common.h:29
MAX
#define MAX(a, b)
Definition utils_common.h:66
points
static struct user_input_data_points ** points
Definition yac.c:158