YetAnotherCoupler 3.2.0_a
Loading...
Searching...
No Matches
test_area.c

A test for different area calculations.

// Copyright (c) 2024 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "test_common.h"
#include "clipping.h"
#include "area.h"
#include "tests.h"
#include "geometry.h"
#define YAC_EARTH_RADIUS (6371.2290)
//#define VERBOSE
static int compare_area(double a, double b);
static void test_area(struct yac_grid_cell cell, double ref_area,
double (*area_func)(struct yac_grid_cell));
int main (void) {
enum yac_edge_type gc_edges[5] =
{YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE};
enum yac_edge_type lat_edges[4] =
{YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
double area;
#ifdef VERBOSE
puts ("----- calculating area -----\n");
#endif
{
/* set the test data over the Equator
0.0 (lon)
0.5 (lat)
/ \
/ \
-0.5 / \ 0.5
0.0 \ / 0.0
\ /
\ /
0.0
-0.5
*/
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){0.5, 0.0, -0.5, 0.0},
(double[4]){0.0, 0.5, 0.0, -0.5}, gc_edges, 4);
area = yac_pole_area ( Quad );
#ifdef VERBOSE
printf ( "pole area of quad over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Quad );
#ifdef VERBOSE
printf ( "huiliers area of quad over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Quad);
}
{
/* Triangle, half of the above Quad */
struct yac_grid_cell Tri =
generate_cell_deg((double[3]){0.5, 0.0, -0.5},
(double[3]){0.0, 0.5, 0.0}, gc_edges, 3);
area = yac_triangle_area ( Tri );
#ifdef VERBOSE
printf ( "ICON area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_pole_area ( Tri );
#ifdef VERBOSE
printf ( "pole area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Tri );
#ifdef VERBOSE
printf ( "huiliers area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Tri);
}
{
/* set the test data over the Equator
-0.5 0.5 (lon)
0.5 ----------- 0.5 (lat)
| |
| |
| |
| |
| |
-0.5 ----------- 0.5
-0.5 -0.5
*/
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){-0.5, 0.5, 0.5, -0.5},
(double[4]){-0.5, -0.5, 0.5, 0.5}, gc_edges, 4);
area = yac_pole_area ( Quad );
#ifdef VERBOSE
printf ( "pole area of quad over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Quad );
#ifdef VERBOSE
printf ( "huiliers area of quad over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
double ref_area = yac_pole_area(Quad);
test_area(Quad, ref_area, yac_pole_area);
test_area(Quad, ref_area, yac_huiliers_area);
yac_free_grid_cell(&Quad);
}
{
/* Half of the above Quad */
struct yac_grid_cell Tri =
generate_cell_deg((double[3]){-0.5, 0.5, 0.5},
(double[3]){-0.5, -0.5, 0.5}, gc_edges, 3);
area = yac_triangle_area ( Tri );
#ifdef VERBOSE
printf ( "ICON area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_pole_area ( Tri );
#ifdef VERBOSE
printf ( "pole area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Tri );
#ifdef VERBOSE
printf ( "huiliers area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Tri);
}
{
struct yac_grid_cell Tri =
generate_cell_deg((double[3]){150.0, 180.0, 180.0},
(double[3]){-60.0, -60.0, -90.0}, gc_edges, 3);
area = yac_triangle_area ( Tri );
#ifdef VERBOSE
printf ( "ICON area of triangle over South Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_pole_area ( Tri );
#ifdef VERBOSE
printf ( "pole area of triangle over South Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Tri );
#ifdef VERBOSE
printf ( "huiliers area of triangle over South Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Tri);
}
{
/* set the test data over the North Pole
90.0 0.0 (lon)
89.5 ----------- 89.5 (lat)
| |
| |
| x |
| |
| |
180.0 ----------- 270.0
89.5 89.5
*/
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){0.0, 90.0, -180.0, -90.0},
(double[4]){89.5, 89.5, 89.5, 89.5}, gc_edges, 4);
area = yac_pole_area ( Quad );
#ifdef VERBOSE
printf ( "pole area of quad over North Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Quad );
#ifdef VERBOSE
printf ( "huiliers area of quad over North Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Quad);
}
{
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){-180.0, -150.0, -150.0, -180.0},
(double[4]){60.0, 60.0, 90.0, 90.0}, gc_edges, 4);
area = yac_pole_area ( Quad );
#ifdef VERBOSE
printf ( "pole area of quad near North Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Quad );
#ifdef VERBOSE
printf ( "huiliers area of quad near North Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Quad);
}
{
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){-180.0, -150.0, -150.0, -180.0},
(double[4]){-90.0, -90.0, -60.0, -60.0}, gc_edges, 4);
#ifdef VERBOSE
printf ( "ICON area of quad near South Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_pole_area ( Quad );
#ifdef VERBOSE
printf ( "pole area of quad near South Pole is %f sqr Km \n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Quad );
#ifdef VERBOSE
printf ( "huiliers area of quad near South Pole is %f sqr Km \n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Quad);
}
{
struct yac_grid_cell Quin =
generate_cell_deg((double[5]){-180.0, -150.0, 0.0, -150.0, -180.0},
(double[5]){60.0, 60.0, 90.0, 90.0, 90.0}, gc_edges, 5);
area = yac_pole_area ( Quin );
#ifdef VERBOSE
printf ( "pole area of quin near North Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Quin );
#ifdef VERBOSE
printf ( "huiliers area of quin near North Pole is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Quin);
}
{
/* Triangle that contains an edge of length zero */
struct yac_grid_cell Tri =
generate_cell_deg((double[3]){-0.5, 0.5, -0.5},
(double[3]){-0.5, -0.5, -0.5}, gc_edges, 3);
area = yac_triangle_area ( Tri );
#ifdef VERBOSE
printf ( "ICON area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_pole_area ( Tri );
#ifdef VERBOSE
printf ( "pole area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area ( Tri );
#ifdef VERBOSE
printf ( "huiliers area of triangle over Equator is %f sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
if ((area < 0.0) || (area > 1e-9))
PUT_ERR("ERROR in yac_huiliers_area for zero size triangle");
yac_free_grid_cell(&Tri);
}
{
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){0, 1, 1, 0},
(double[4]){0, 0, 1, 1}, gc_edges, 4);
struct yac_grid_cell Tri =
generate_cell_deg((double[3]){1, 1, 1},
(double[3]){0, 0, 1}, gc_edges, 3);
double base_area = yac_huiliers_area(Quad);
yac_free_grid_cell(&Quad);
yac_free_grid_cell(&Tri);
/* A dx of 1e-8 deg corresponds to a dx of ~1 mm on the Earth surface */
double dx[] = {1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9};
for (size_t i = 0; i < sizeof(dx) / sizeof(dx[0]); ++i) {
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){0, 1.0-dx[i], 1, 0},
(double[4]){0, 0, 1, 1}, gc_edges, 4);
struct yac_grid_cell Tri =
generate_cell_deg((double[3]){1.0-dx[i], 1, 1},
(double[3]){0, 0, 1}, gc_edges, 3);
double partial_area = yac_huiliers_area(Quad);
partial_area += yac_huiliers_area(Tri);
yac_free_grid_cell(&Quad);
yac_free_grid_cell(&Tri);
if (compare_area(base_area, partial_area))
PUT_ERR("yac_huiliers_area computed wrong area\n");
#ifdef VERBOSE
printf ("accuracy check for (base_area - partial_area) %8.1e sqr km for dx %.1e\n",
base_area - partial_area, dx[i]);
#endif
}
}
{
struct yac_grid_cell Tri =
generate_cell_deg(
(double[3]){0.09817477/YAC_RAD, 0.09817477/YAC_RAD, 0.098174496/YAC_RAD},
(double[3]){0.353325247/YAC_RAD, 0.353335501/YAC_RAD, 0.353335609/YAC_RAD},
gc_edges, 3);
area = yac_pole_area(Tri);
if (area < 0)
PUT_ERR("pole_area computed wrong area\n");
#ifdef VERBOSE
printf ("pole accuracy check of very small area %.2e sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
area = yac_huiliers_area(Tri);
if (area < 0)
PUT_ERR("huiliers_area computed wrong area\n");
#ifdef VERBOSE
printf ("huiliers accuracy check of very small area %.2e sqr km\n", area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS );
#endif
yac_free_grid_cell(&Tri);
}
{
double edge_length_deg = 1.0e-8; // ~1.11 mm
#ifdef VERBOSE
printf ("\n accuracy check for a unit square\n");
printf (" edge length [m] | plain | simple | pole | huiliers\n");
printf (" ----------------+--------------------+--------------------+--------------------+--------------------\n");
#endif
for ( unsigned i = 1; i < 11; ++i ) {
double edge_length_m =
edge_length_deg * YAC_RAD * YAC_EARTH_RADIUS * 1.0e3;
struct yac_grid_cell Quad =
generate_cell_deg((double[4]){0, edge_length_deg, edge_length_deg, 0},
(double[4]){0, 0, edge_length_deg, edge_length_deg},
gc_edges, 4);
double t_area = edge_length_m * edge_length_m;
double p_area = yac_planar_3dcell_area(Quad) * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS * 1.0e6;
double s_area = yac_pole_area(Quad) * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS * 1.0e6;
double h_area = yac_huiliers_area(Quad) * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS * 1.0e6;
#ifndef VERBOSE
if (0)
#endif
printf (" %.8e m|%.8e sqr m|%.8e sqr m|%.8e sqr m|%.8e sqr m\n",
edge_length_m, p_area, t_area, s_area, h_area);
edge_length_deg *= 10.0;
yac_free_grid_cell(&Quad);
}
}
{ // cell with two edges (one great circle and one circle of latitude)
for (int y = -90; y <= 90; y += 5) {
struct yac_grid_cell cell =
generate_cell_deg((double[2]){0.0, 20.0},
(double[2]){(double)y, (double)y},
(enum yac_edge_type[]){
YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE}, 2);
// compute reference area
double ref_area;
{
double coordinates_x[3] = {0.0, 20.0, 0.0};
double coordinates_y[3] = {(double)y, (double)y, 0};
if ((y == -90) || (y == 90) || (y == 0)) {
ref_area = 0.0;
} else {
if (y < 0) {
ref_area = (coordinates_x[1]*YAC_RAD - coordinates_x[0]*YAC_RAD) *
(1 - cos(M_PI_2 + coordinates_y[0]*YAC_RAD));
coordinates_y[2] = -90.0;
} else {
ref_area = (coordinates_x[1]*YAC_RAD - coordinates_x[0]*YAC_RAD) *
(1 - cos(M_PI_2 - coordinates_y[0]*YAC_RAD));
coordinates_y[2] = 90.0;
}
struct yac_grid_cell tmp_cell =
generate_cell_deg(coordinates_x, coordinates_y, gc_edges, 3);
// ref_area *= YAC_EARTH_RADIUS * YAC_EARTH_RADIUS;
ref_area -= yac_huiliers_area(tmp_cell);
yac_free_grid_cell(&tmp_cell);
}
}
test_area(cell, ref_area, yac_pole_area);
test_area(cell, ref_area, yac_huiliers_area);
yac_free_grid_cell(&cell);
}
}
{ // cell with lon lat edges and one edge on the equator
double ref_area = sin(0.5 * YAC_RAD) * (0.5 * YAC_RAD);
// ref_area *= YAC_EARTH_RADIUS * YAC_EARTH_RADIUS;
double coordinates_x[2][4] = {{-0.5, 0.0, 0.0, -0.5},
{ 0.0, 0.5, 0.5, 0.0}};
double coordinates_y[2][4] = {{-0.5, -0.5, 0.0, 0.0},
{ 0.0, 0.0, 0.5, 0.5}};
for (int x = 0; x < 2; ++x) {
for (int y = 0; y < 2; ++y) {
struct yac_grid_cell cell =
generate_cell_deg(coordinates_x[x], coordinates_y[y], lat_edges, 4);
test_area(cell, ref_area, yac_pole_area);
test_area(cell, ref_area, yac_huiliers_area);
yac_free_grid_cell(&cell);
}
}
}
{ // cell with lon lat edges (both lat edges on the different side of the
// equator)
double ref_area = (2.0 * sin(0.5 * YAC_RAD)) * (1.0 * YAC_RAD);
// ref_area *= YAC_EARTH_RADIUS * YAC_EARTH_RADIUS;
struct yac_grid_cell cell =
generate_cell_deg((double[4]){-0.5, 0.5, 0.5, -0.5},
(double[4]){-0.5, -0.5, 0.5, 0.5}, lat_edges, 4);
test_area(cell, ref_area, yac_pole_area);
test_area(cell, ref_area, yac_huiliers_area);
yac_free_grid_cell(&cell);
}
{ // cell with lon lat edges (both lat edges on the same side of the equator)
double ref_area = sin(1.0 * YAC_RAD);
ref_area -= sin(0.5 * YAC_RAD);
ref_area *= 1.0 * YAC_RAD;
// ref_area *= YAC_EARTH_RADIUS * YAC_EARTH_RADIUS;
for (int i = 0; i < 2; ++i) {
double coordinates_x[4] = {-0.5, 0.5, 0.5, -0.5};
double coordinates_y[2][4] = {{ 0.5, 0.5, 1.0, 1.0},
{-0.5,-0.5,-1.0, -1.0}};
enum yac_edge_type edge_types[4] =
{YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
struct yac_grid_cell cell =
generate_cell_deg(coordinates_x, coordinates_y[i], edge_types, 4);
test_area(cell, ref_area, yac_pole_area);
test_area(cell, ref_area, yac_huiliers_area);
yac_free_grid_cell(&cell);
}
}
{ // concave cell
double intersections[2][3];
if (!intersect(YAC_GREAT_CIRCLE_EDGE, 20.0, 59.0, 0.0, 80.0,
YAC_LAT_CIRCLE_EDGE, 0.0, 60.0, 20.0, 60.0,
intersections[0]))
return EXIT_FAILURE;
if (!intersect(YAC_GREAT_CIRCLE_EDGE, 20.0, 59.0, -20.0, 59.0,
YAC_LAT_CIRCLE_EDGE, 0.0, 60.0, 20.0, 60.0,
intersections[1]))
return EXIT_FAILURE;
double intersections_lon[2], intersections_lat[2];
XYZtoLL(intersections[0], &(intersections_lon[0]), &(intersections_lat[0]));
XYZtoLL(intersections[1], &(intersections_lon[1]), &(intersections_lat[1]));
// this cell is convex
struct yac_grid_cell partial_ref_cell =
generate_cell_deg(
(double[]){20.0, intersections_lon[0]/YAC_RAD,
intersections_lon[1]/YAC_RAD},
(double[]){59.0, 60.0, 60.0},
(enum yac_edge_type[]){YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 3);
double partial_area = yac_pole_area(partial_ref_cell);
if (compare_area(partial_area,
yac_huiliers_area(partial_ref_cell)))
PUT_ERR("pole_area != huiliers_area for convex cell");
double ref_area = 2.0 * partial_area;
struct yac_grid_cell cells[2] =
{generate_cell_deg(
(double[]){20.0,
intersections_lon[0]/YAC_RAD,
intersections_lon[1]/YAC_RAD,
-intersections_lon[1]/YAC_RAD,
-intersections_lon[0]/YAC_RAD, -20.0},
(double[]){59.0, 60.0, 60.0, 60.0, 60.0, 59.0},
(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_GREAT_CIRCLE_EDGE}, 6),
generate_cell_deg(
(double[]){20.0,
intersections_lon[0]/YAC_RAD,
-intersections_lon[0]/YAC_RAD,
-20.0,
-intersections_lon[1]/YAC_RAD,
intersections_lon[1]/YAC_RAD},
(double[]){59.0, 60.0, 60.0, 59.0, 60.0, 60.0},
(enum yac_edge_type[]){
YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}, 6)};
for (int i = 0; i < 2; ++i) {
test_area(cells[i], ref_area, yac_pole_area);
test_area(cells[i], ref_area, yac_huiliers_area);
}
yac_free_grid_cell(&(cells[0]));
yac_free_grid_cell(&(cells[1]));
yac_free_grid_cell(&partial_ref_cell);
}
{ // case extracted from ICON example, which requires a high accuracy in the
// area computation
struct yac_grid_cell cells[] = {
generate_cell_deg((double[]){-151.1489316535362661,
-151.7087167856307133,
-150.6174180645230933},
(double[]){ 10.6333746562670566,
11.5836995038666437,
11.6243355661022782}, gc_edges, 3),
generate_cell_deg((double[]){-152.8132296086113513,
-153.8786613849353557,
-151.6929566531634066},
(double[]){ 13.4512385235805798,
11.4650630285209161,
11.5568484830228009}, gc_edges, 3),
generate_cell_deg((double[]){-152.8132296086113513,
-151.6929566531634066,
-150.6265735326831532},
(double[]){ 13.4512385235805798,
11.5568484830228009,
13.5446528932721275}, gc_edges, 3),
generate_cell_deg((double[]){-153.8786613849353557,
-151.6929566531634066,
-152.7627337461449031},
(double[]){ 11.4650630285209161,
11.5568484830228009,
9.5827293336264265}, gc_edges, 3),
generate_cell_deg((double[]){-151.6929566531634066,
-152.7627337461449031,
-150.5840076321059371},
(double[]){ 11.5568484830228009,
9.5827293336264265,
9.6520960946755832}, gc_edges, 3),
generate_cell_deg((double[]){-151.6929566531634066,
-150.6265735326831532,
-149.5139953075027108},
(double[]){ 11.5568484830228009,
13.5446528932721275,
11.6284420125700851}, gc_edges, 3),
generate_cell_deg((double[]){-151.6929566531634066,
-149.5139953075027108,
-150.5840076321059371},
(double[]){ 11.5568484830228009,
11.6284420125700851,
9.6520960946755832}, gc_edges, 3)};
size_t num_cells = sizeof(cells) / sizeof(cells[0]);
struct yac_grid_cell overlap_cells[sizeof(cells) / sizeof(cells[0]) - 1];
for (size_t i = 0; i < num_cells-1; ++i)
yac_init_grid_cell(&overlap_cells[i]);
yac_cell_clipping(num_cells - 1, &cells[1], cells[0], &overlap_cells[0]);
double tgt_area = yac_huiliers_area(cells[0]);
double area_tol = tgt_area * 1e-8;
double diff_area = tgt_area;
for (size_t i = 0; i < num_cells-1; ++i)
diff_area -= yac_huiliers_area(overlap_cells[i]);
if (fabs(diff_area) > area_tol)
PUT_ERR("yac_huiliers_area is too inaccurate");
for (size_t i = 0; i < num_cells; ++i)
yac_free_grid_cell(&cells[i]);
for (size_t i = 0; i < num_cells-1; ++i)
yac_free_grid_cell(&overlap_cells[i]);
}
{ // check for a bug found by Uwe Schulzweida
struct yac_grid_cell Tri =
generate_cell_deg((double[3]){0.5, 0.5, -0.5},
(double[3]){0.0, 0.0, 0.0}, gc_edges, 3);
if (yac_triangle_area(Tri) != 0.0)
PUT_ERR("ERROR in yac_triangle_area");
yac_free_grid_cell(&Tri);
Tri = generate_cell_deg((double[3]){0.5, -0.5, 0.5},
(double[3]){0.0, 0.0, 0.0}, gc_edges, 3);
if (yac_triangle_area(Tri) != 0.0)
PUT_ERR("ERROR in yac_triangle_area");
yac_free_grid_cell(&Tri);
Tri = generate_cell_deg((double[3]){-0.5, 0.5, 0.5},
(double[3]){0.0, 0.0, 0.0}, gc_edges, 3);
if (yac_triangle_area(Tri) != 0.0)
PUT_ERR("ERROR in yac_triangle_area");
yac_free_grid_cell(&Tri);
}
return TEST_EXIT_CODE;
}
static void test_area(struct yac_grid_cell cell, double ref_area,
double (*area_func)(struct yac_grid_cell)) {
double mem_dummy_[128][3];
enum yac_edge_type edge_dummy[128];
struct yac_grid_cell test_cell = {.coordinates_xyz = mem_dummy_,
.edge_type = edge_dummy,
.num_corners = cell.num_corners};
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];
}
double area = area_func(test_cell);
if (compare_area(area, ref_area)) PUT_ERR("ERROR: wrong area\n");
}
}
}
static int compare_area(double a, double b) {
double diff = fabs(a - b);
double tol = MAX(MIN(a,b)*1e-6, YAC_AREA_TOL);
if (diff < tol) return 0;
else return (a > b) - (a < b);
}
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:396
yac_pole_area
double yac_pole_area(struct yac_grid_cell cell)
Calculate the area of a cell in a 3d plane on a unit sphere.
Definition area.c:279
yac_planar_3dcell_area
double yac_planar_3dcell_area(struct yac_grid_cell cell)
Area calculation on a unit sphere of a planar polygon in 3D.
Definition area.c:364
yac_triangle_area
double yac_triangle_area(struct yac_grid_cell cell)
Calculate the area of a triangle on a unit sphere.
Definition area.c:21
area.h
Structs and interfaces for area calculations.
YAC_AREA_TOL
#define YAC_AREA_TOL
Definition area.h:23
tol
static double const tol
Definition check_overlap.c:11
yac_cell_clipping
void yac_cell_clipping(size_t N, struct yac_grid_cell *source_cell, struct yac_grid_cell target_cell, struct yac_grid_cell *overlap_buffer)
cell clipping to get the cells describing the intersections
Definition clipping.c:1078
YAC_RAD
#define YAC_RAD
Definition geometry.h:30
XYZtoLL
static void XYZtoLL(double const p_in[], double *lon, double *lat)
Definition geometry.h:318
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
yac_replay.cells
cells
Definition yac_replay.py:117
yac_replay.y
y
Definition yac_replay.py:57
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
MIN
#define MIN(a, b)
Definition utils_common.h:70
MAX
#define MAX(a, b)
Definition utils_common.h:66
main
int main(int argc, char **argv)
Definition weights2vtk.c:147
Generated on Mon Apr 8 2024 08:40:58 for YetAnotherCoupler by doxygen 1.9.8