test_lat_clipping.c

Go to the documentation of this file.

// Copyright (c) 2024 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
#include "grids/basic_grid.h"
#include "clipping.h"
#include "geometry.h"
#include "tests.h"
#include "area.h"
#include "test_common.h"
 
static void utest_lat_clipping(struct yac_grid_cell cells,
                              double deg_lat_bounds[2],
                              struct yac_grid_cell * ref_cells,
                              unsigned num_ref_cells);
 
static int utest_compare_cells(struct yac_grid_cell a, struct yac_grid_cell b);
 
static void utest_lat_clipping_degenerate_cells(void);
 
static double utest_get_intersection(
  enum yac_edge_type edge_type, double points[2][2], double lat) {
 
  double intersection[3], intersection_lon, intersection_lat;
 
  if (!utest_intersect(edge_type, points[0][0], points[0][1],
                            points[1][0], points[1][1],
                 YAC_LAT_CIRCLE_EDGE, points[0][0], lat,
                                  points[1][0], lat, intersection))
    exit(EXIT_FAILURE);
  XYZtoLL(intersection, &intersection_lon, &intersection_lat);
  return intersection_lon / YAC_RAD;
}
 
int main (void) {
 
  enum yac_edge_type gc_edges[] = {
    YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
    YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
    YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
    YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE};
  enum yac_edge_type latlon_edges[] = {
    YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
  // edges of a lon-lat cell whose lat edge collapsed at a pole
  // (the pole is the first corner)
  enum yac_edge_type lonlat_pole_edges[] = {
    YAC_LON_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
  // edges of a cell whose boundary is completely on a circle of latitude
  enum yac_edge_type lat3_edges[] = {
    YAC_LAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE};
 
  double temp_lon[16];
 
  // a gaussian cell outside the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){-5,-5,5,5}, latlon_edges, 4),
    (double[2]){10, 20}, NULL, 0);
 
  // a gaussian cell outside the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){-5,-5,10,10}, latlon_edges, 4),
    (double[2]){10, 20}, NULL, 0);
 
  // a gaussian cell outside the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){20,20,30,30}, latlon_edges, 4),
    (double[2]){10, 20}, NULL, 0);
 
  // a gaussian cell outside the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){80,80,90,90}, latlon_edges, 4),
    (double[2]){10, 20}, NULL, 0);
 
  // a gaussian cell outside the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){-80,-80,-90,-90}, latlon_edges, 4),
    (double[2]){10, 20}, NULL, 0);
 
  // a gaussian cell outside the lat bound (one bound is at a pole)
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){80, 90}, NULL, 0);
 
  // a gaussian cell outside the lat bound (both bounds are identical)
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){0, 0}, NULL, 0);
 
  // a gaussian cell outside the lat bound (both bounds are identical)
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){-10,-10,-20,-20}, latlon_edges, 4),
    (double[2]){0, 0}, NULL, 0);
 
  // a gaussian cell outside the lat bound (both bounds are identical)
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){-5,-5,5,5}, latlon_edges, 4),
    (double[2]){0, 0}, NULL, 0);
 
  // a gaussian cell overlapping with the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){-10,-10,10,10}, latlon_edges, 4),
    (double[2]){-5, 5},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){-5,5,5,-5},
       (double[4]){-5,-5,5,5}, latlon_edges, 4)}, 1);
 
  // a gaussian cell overlapping with the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){12, 15},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){12,12,15,15}, latlon_edges, 4)}, 1);
 
  // a gaussian cell overlapping with the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){15, 25},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){15,15,20,20}, latlon_edges, 4)}, 1);
 
  // a gaussian cell overlapping with the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){15, 5},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){15,15,10,10}, latlon_edges, 4)}, 1);
 
  // a gaussian cell overlapping with the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){5, 25},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){10,10,20,20}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){90, 80}, NULL, 0);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){90,90,70,70}, latlon_edges, 4),
    (double[2]){90, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){90,90,80,80}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){90,90,85,85}, latlon_edges, 4),
    (double[2]){90, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){90,90,85,85}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){88,88,83,83}, latlon_edges, 4),
    (double[2]){90, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){88,88,83,83}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){88,88,83,83}, latlon_edges, 4),
    (double[2]){90, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){88,88,83,83}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){88,88,70,70}, latlon_edges, 4),
    (double[2]){90, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){88,88,80,80}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){10,10,20,20}, latlon_edges, 4),
    (double[2]){-90, -80}, NULL, 0);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-90,-90,-70,-70}, latlon_edges, 4),
    (double[2]){-90, -80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-90,-90,-80,-80}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-90,-90,-85,-85}, latlon_edges, 4),
    (double[2]){-90, -80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-90,-90,-85,-85}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-88,-88,-83,-83}, latlon_edges, 4),
    (double[2]){-90, -80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-88,-88,-83,-83}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-88,-88,-83,-83}, latlon_edges, 4),
    (double[2]){-90, -80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-88,-88,-83,-83}, latlon_edges, 4)}, 1);
 
  // a gaussian cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-88,-88,-70,-70}, latlon_edges, 4),
    (double[2]){-90, -80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[4]){40,45,45,40},
       (double[4]){-88,-88,-80,-80}, latlon_edges, 4)}, 1);
 
  // a gaussian cell that touches the north pole and partially overlaps
  // with the bounds
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[4]){0,5,5,0},
       (double[4]){90,90,85,85}, latlon_edges, 4),
    (double[2]){89.9, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
        (double[4]){0,5,5,0},
        (double[4]){89.9,89.9,85,85}, latlon_edges, 4)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[3]){10,20,15},
       (double[3]){30,30,40}, gc_edges, 3),
    (double[2]){35, 50},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){15,
                  utest_get_intersection(
                     YAC_GREAT_CIRCLE_EDGE,
                     (double[2][2]){{15, 40}, {10, 30}}, 35),
                  utest_get_intersection(
                     YAC_GREAT_CIRCLE_EDGE,
                     (double[2][2]){{15, 40}, {20, 30}}, 35)},
       (double[]){40,35,35},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 3)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[3]){10,20,15},
       (double[3]){30,30,40}, gc_edges, 3),
    (double[2]){35, 25},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){20, 10,
                  utest_get_intersection(
                     YAC_GREAT_CIRCLE_EDGE,
                     (double[2][2]){{15, 40}, {10, 30}}, 35),
                  utest_get_intersection(
                     YAC_GREAT_CIRCLE_EDGE,
                     (double[2][2]){{15, 40}, {20, 30}}, 35)},
       (double[]){30,30,35,35},
       (enum yac_edge_type[]){YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
                              YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 4)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){30,30,40}, gc_edges, 3),
    (double[2]){45, 25},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){30,30,40}, gc_edges, 3)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[3]){10,20,15},
       (double[3]){30,30,40}, gc_edges, 3),
    (double[2]){32, 38},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{15, 40}, {20, 30}}, 32),
                  utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{15, 40}, {20, 30}}, 38),
                  utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{15, 40}, {10, 30}}, 38),
                  utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{15, 40}, {10, 30}}, 32)},
       (double[]){32,38,38,32},
       (enum yac_edge_type[]){YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE,
                              YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE}, 4)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[3]){-20,20,0},
       (double[3]){80,80,85}, gc_edges, 3),
    (double[2]){80.001, 89},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,
                  utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{0, 85}, {-20, 80}}, 80.001),
                  -fabs(utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{-20, 80}, {20, 80}}, 80.001)),
                  fabs(utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{-20, 80}, {20, 80}}, 80.001)),
                  utest_get_intersection(
                    YAC_GREAT_CIRCLE_EDGE,
                    (double[2][2]){{0, 85}, {20, 80}}, 80.001),},
       (double[]){85,80.001,80.001,80.001,80.001},
       (enum yac_edge_type[]){YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE,
                              YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE,
                              YAC_GREAT_CIRCLE_EDGE}, 5)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){82,82,88}, gc_edges, 3),
    (double[2]){90, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){82,82,88}, gc_edges, 3)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){-82,-82,-88}, gc_edges, 3),
    (double[2]){90, 80}, NULL, 0);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){-82,-82,-88}, gc_edges, 3),
    (double[2]){-90, -80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){-82,-82,-88}, gc_edges, 3)}, 1);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){82,82,88}, gc_edges, 3),
    (double[2]){-90, -80}, NULL, 0);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){82,82,88}, gc_edges, 3),
    (double[2]){90, 90}, NULL, 0);
 
  // a triangle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){10,20,15},
       (double[]){82,82,88}, gc_edges, 3),
    (double[2]){-90, -90}, NULL, 0);
 
  // a triangle cell covering the pole
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){80,80,80}, gc_edges, 3),
    (double[2]){90, 70},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){80,80,80}, gc_edges, 3)}, 1);
 
  // a triangle cell covering the pole
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){80,80,80}, gc_edges, 3),
    (double[2]){60, 70}, NULL, 0);
 
  // a triangle cell covering the pole
  temp_lon[0] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
  temp_lon[1] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){85,70,70}, gc_edges, 3),
    (double[2]){90, 80},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,temp_lon[0],temp_lon[1]},
       (double[]){85,80,80},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 3)}, 1);
 
  // a triangle cell covering the pole
  temp_lon[0] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 72);
  temp_lon[1] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{120, 70}, {240, 70}}, 72);
  temp_lon[1] = (temp_lon[1] < 0)?(temp_lon[1] + 360):(temp_lon[1]);
  temp_lon[1] = (temp_lon[1] < 180)?(temp_lon[1]):(360 - temp_lon[1]);
  temp_lon[2] = 360 - temp_lon[1];
  temp_lon[3] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 72);
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){85,70,70}, gc_edges, 3),
    (double[2]){90, 72},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,temp_lon[0],temp_lon[1],temp_lon[2],temp_lon[3]},
       (double[]){85,72,72,72,72},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
         YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 5)}, 1);
 
  // a triangle cell covering the pole
  temp_lon[0] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
  temp_lon[1] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){85,70,70}, gc_edges, 3),
    (double[2]){80, 65},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){temp_lon[0],120,240,temp_lon[1]},
       (double[]){80,70,70,80},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
         YAC_LAT_CIRCLE_EDGE}, 4)}, 1);
 
  // a triangle cell covering the pole
  temp_lon[0] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 84);
  temp_lon[1] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
  temp_lon[2] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
  temp_lon[3] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 84);
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){85,70,70}, gc_edges, 3),
    (double[2]){80, 84},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){temp_lon[0],temp_lon[1],temp_lon[2],temp_lon[3]},
       (double[]){84,80,80,84},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
         YAC_LAT_CIRCLE_EDGE}, 4)}, 1);
 
  // a triangle cell covering the pole
  temp_lon[0] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 80);
  temp_lon[1] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {120, 70}}, 72);
  temp_lon[2] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{240, 70}, {120, 70}}, 72);
  temp_lon[2] = (temp_lon[2] < 0)?(temp_lon[2] + 360):(temp_lon[2]);
  temp_lon[2] = (temp_lon[2] < 180)?(temp_lon[2]):(360 - temp_lon[2]);
  temp_lon[3] = 360 - temp_lon[2];
  temp_lon[4] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 72);
  temp_lon[5] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 85}, {240, 70}}, 80);
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){85,70,70}, gc_edges, 3),
    (double[2]){80, 72},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){temp_lon[0],temp_lon[1],temp_lon[2],
                  temp_lon[3],temp_lon[4],temp_lon[5]},
       (double[]){80,72,72,72,72,80},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
         YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE}, 6)}, 1);
 
  // a triangle cell touching the pole
  temp_lon[0] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 90}, {120, 80}}, 85);
  temp_lon[1] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 90}, {240, 80}}, 85);
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,120,240},
       (double[]){90,80,80}, gc_edges, 3),
    (double[2]){90, 85},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,temp_lon[0],temp_lon[1]},
       (double[]){90,85,85},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 3)}, 1);
 
  // a triangle cell touching the pole
  temp_lon[0] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{0, 0}, {10, 5}}, 2.5);
  temp_lon[1] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{10, -5}, {10, 5}}, 2.5);
  temp_lon[2] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{10, -5}, {10, 5}}, -2.5);
  temp_lon[3] =
    utest_get_intersection(YAC_GREAT_CIRCLE_EDGE, (double[2][2]){{10, -5}, {0, 0}}, -2.5);
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,10,10},
       (double[]){0,5,-5}, gc_edges, 3),
    (double[2]){2.5, -2.5},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,temp_lon[0],temp_lon[1],temp_lon[2],temp_lon[3]},
       (double[]){0,2.5,2.5,-2.5,-2.5},
       (enum yac_edge_type[]){
         YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
         YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 5)}, 1);
 
  // A lon-lat cell whose upper lat edge collapsed into a single point at the
  // pole. Such a cell is stored with three corners (the pole being the first
  // one here), but is still classified as a lon-lat cell.
 
  // the cell is completely within the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,0,90},
       (double[]){90,70,70}, lonlat_pole_edges, 3),
    (double[2]){60, 90},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,0,90},
       (double[]){90,70,70}, lonlat_pole_edges, 3)}, 1);
 
  // the cell is outside of the lat bound
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,0,90},
       (double[]){90,70,70}, lonlat_pole_edges, 3),
    (double[2]){30, 60}, NULL, 0);
 
  // only the lower lat bound overlaps with the cell
  // (the pole remains part of the overlap)
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,0,90},
       (double[]){90,70,70}, lonlat_pole_edges, 3),
    (double[2]){75, 90},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,0,90},
       (double[]){90,75,75}, lonlat_pole_edges, 3)}, 1);
 
  // both lat bounds overlap with the cell
  // (the pole is clipped off, which turns the triangle into a lon-lat cell
  //  with four corners)
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,0,90},
       (double[]){90,70,70}, lonlat_pole_edges, 3),
    (double[2]){75, 85},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,90,90,0},
       (double[]){75,75,85,85}, latlon_edges, 4)}, 1);
 
  // A cell covering the pole, whose boundary consists of three edges on the
  // same circle of latitude (covering 170, 170, and 20 degree of longitude).
  // For two consecutive of these edges the longitude difference between their
  // first and their last vertex is only 20 degree, even though the two edges
  // together cover 340 degree.
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,170,340},
       (double[]){45,45,45}, lat3_edges, 3),
    (double[2]){0, 90},
    (struct yac_grid_cell[]){
      utest_generate_cell_deg(
       (double[]){0,170,340},
       (double[]){45,45,45}, lat3_edges, 3)}, 1);
 
  // empty great circle cell
  utest_lat_clipping(
    utest_generate_cell_deg(
       (double[]){0,0,0},
       (double[]){-5,0,5}, gc_edges, 3),
    (double[2]){2.5, -2.5}, NULL, 0);
 
  // cells with too few corners
  utest_lat_clipping_degenerate_cells();
 
  return TEST_EXIT_CODE;
}
 
static void utest_lat_clipping_degenerate_cells(void) {
 
  enum yac_edge_type latlon_edges[4] = {
    YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE,
    YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
  double rad_lat_bounds[2] = {-10.0 * YAC_RAD, 10.0 * YAC_RAD};
 
  struct yac_grid_cell overlap_buffer[2];
  for (int i = 0; i < 2; ++i) yac_init_grid_cell(overlap_buffer + i);
 
  // fill the overlap buffer with non-empty overlaps
  struct yac_grid_cell cells[2] = {
    utest_generate_cell_deg(
      (double[4]){-5,5,5,-5}, (double[4]){-5,-5,5,5}, latlon_edges, 4),
    utest_generate_cell_deg(
      (double[4]){5,15,15,5}, (double[4]){-5,-5,5,5}, latlon_edges, 4)};
  yac_cell_lat_clipping(2, cells, rad_lat_bounds, overlap_buffer);
  for (int i = 0; i < 2; ++i) {
    if (overlap_buffer[i].num_corners == 0)
      PUT_ERR("ERROR: test setup failed (expected a non-empty overlap)\n");
    yac_free_grid_cell(cells + i);
  }
 
  // clip degenerate cells using the same (non-empty) overlap buffer
  double dummy_coordinates_xyz[1][3];
  LLtoXYZ_deg(0.0, 0.0, dummy_coordinates_xyz[0]);
  enum yac_edge_type dummy_edge_type[1] = {YAC_GREAT_CIRCLE_EDGE};
  struct yac_grid_cell degenerate_cells[2];
  for (int i = 0; i < 2; ++i)
    degenerate_cells[i] =
      (struct yac_grid_cell){.coordinates_xyz = dummy_coordinates_xyz,
                             .edge_type       = dummy_edge_type,
                             .num_corners     = (size_t)i,
                             .array_size      = 1};
  yac_cell_lat_clipping(2, degenerate_cells, rad_lat_bounds, overlap_buffer);
  for (int i = 0; i < 2; ++i)
    if (overlap_buffer[i].num_corners != 0)
      PUT_ERR("ERROR: cell with too few corners does not result in an "
              "empty overlap\n");
 
  for (int i = 0; i < 2; ++i) yac_free_grid_cell(overlap_buffer + i);
}
 
static int utest_compare_cells(struct yac_grid_cell a, struct yac_grid_cell b) {
 
  // Trigonometric functions can be very inaccurate depending on the compiler
  // and the optimsation level. Therefore, we use angle_tol instead of
  // yac_angle_tol in this routine.
  double const angle_tol = 1e-6;
 
  if ((yac_grid_cell_area(a) < angle_tol * angle_tol) &&
      (yac_grid_cell_area(b) < angle_tol * angle_tol)) return 0;
 
  if (a.num_corners != b.num_corners) return 1;
 
  if (a.num_corners == 0) return 0;
 
  for (int order = -1; order <= 1; order += 2) {
    for (int start = 0; start < (int)(a.num_corners); ++start) {
 
      int differences = 0;
 
      for (int i = 0; i < (int)(a.num_corners); ++i) {
 
        int j =
          ((int)(a.num_corners) + start + order * i) % (int)(a.num_corners);
 
        if (get_vector_angle(
              a.coordinates_xyz[i], b.coordinates_xyz[j]) > angle_tol)
          ++differences;
 
        else if ((a.edge_type[i] == YAC_LAT_CIRCLE_EDGE) !=
                 (b.edge_type[
                    (((int)(a.num_corners))+j-(order<0))%
                    ((int)(a.num_corners))] == YAC_LAT_CIRCLE_EDGE))
          ++differences;
      }
 
      if (!differences) return 0;
    }
  }
 
  return 1;
}
 
static void utest_lat_clipping(struct yac_grid_cell cell,
                              double deg_lat_bounds[2],
                              struct yac_grid_cell * ref_cells,
                              unsigned num_ref_cells) {
 
  double const area_tol = 1e-8 * yac_grid_cell_area(cell);
 
  double rad_lat_bounds[2] =
    {deg_lat_bounds[0] * YAC_RAD, deg_lat_bounds[1] * YAC_RAD};
 
  // the arrays are initialised in order to get a reproducible behaviour in
  // case the clipping accesses more corners than the cell actually has
  double mem_dummy_[128][3] = {{0.0}};
  enum yac_edge_type edge_dummy[128] = {YAC_GREAT_CIRCLE_EDGE};
  struct yac_grid_cell overlap_cell,
                   test_cell =
        (struct yac_grid_cell){.coordinates_xyz = mem_dummy_,
                           .edge_type = edge_dummy,
                           .num_corners = cell.num_corners};
 
  yac_init_grid_cell(&overlap_cell);
 
  for (int pole = 0; pole < 2; ++pole) {
    for (int order = -1; order <= 1; order += 2) {
      for (int start = 0; start < (int)(cell.num_corners); ++start) {
        for (int i = 0; i < ((int)(cell.num_corners)); ++i) {
          int j = (((int)(cell.num_corners))+start+i*order)%
                  ((int)(cell.num_corners));
          test_cell.coordinates_xyz[i][0] = cell.coordinates_xyz[j][0];
          test_cell.coordinates_xyz[i][1] = cell.coordinates_xyz[j][1];
          test_cell.coordinates_xyz[i][2] = cell.coordinates_xyz[j][2];
          j = (((int)(cell.num_corners)) + j - (order < 0))%
              ((int)(cell.num_corners));
          test_cell.edge_type[i] = cell.edge_type[j];
        }
 
        for (int bound_tol_a = -1; bound_tol_a <= 1; ++bound_tol_a) {
          for (int bound_tol_b = -1; bound_tol_b <= 1; ++bound_tol_b) {
            for (int bound_order = 0; bound_order <= 1; ++bound_order) {
              double test_rad_lat_bounds[2] = {
                MIN(M_PI_2,
                  MAX(-M_PI_2,
                      rad_lat_bounds[bound_order] +
                      (double)bound_tol_a * yac_angle_low_tol)),
                MIN(M_PI_2,
                  MAX(-M_PI_2,
                    rad_lat_bounds[bound_order^1] +
                    (double)bound_tol_b * yac_angle_low_tol))};
              yac_cell_lat_clipping(
                1, &test_cell,test_rad_lat_bounds, &overlap_cell);
              if (num_ref_cells == 0) {
                if ((overlap_cell.num_corners != 0) &&
                    (yac_grid_cell_area(overlap_cell) > area_tol))
                  PUT_ERR("ERROR: wrong clipping cell\n");
              } else {
                int match = 0;
                for (int i = 0; i < (int)num_ref_cells; ++i)
                  match |= !utest_compare_cells(overlap_cell, ref_cells[i]);
                if (!match)
                  PUT_ERR("ERROR: wrong clipping cell\n");
              }
            }
          }
        }
      }
    }
 
    // switch to the other pole
    for (int i = 0; i < (int)(cell.num_corners); ++i)
      cell.coordinates_xyz[i][2] *= -1.0;
    rad_lat_bounds[0] *= -1.0, rad_lat_bounds[1] *= -1.0;
    for (unsigned i = 0; i < num_ref_cells; ++i)
      for (int j = 0; j < (int)(ref_cells[i].num_corners); ++j)
        ref_cells[i].coordinates_xyz[j][2] *= -1.0;
  }
 
  yac_free_grid_cell(&cell);
  for (unsigned i = 0; i < num_ref_cells; ++i)
    yac_free_grid_cell(&ref_cells[i]);
  yac_free_grid_cell(&overlap_cell);
}