clipping update

config/default

@@ -28,8 +28,6 @@ case "${HOSTNAME}" in
                     --with-udunits2=$HOME/local/udunits-2.1.24 \
                     --with-proj=/opt/local \
                     --with-curl=/opt/local \
-                    --with-libxml2=/usr \
-                    --with-magics=/Users/m214003/local/magics-2.14.9 \
 	            CC=icc CFLAGS="-g -Wall -O3 -march=native -openmp -fp-model precise" \
                     LIBS="-L/opt/local/lib -lopenjpeg"
 	;;
@@ -45,8 +43,6 @@ case "${HOSTNAME}" in
                     --with-udunits2=$HOME/local/udunits-2.1.24 \
                     --with-proj=/opt/local \
                     --with-curl=/opt/local \
-                    --with-libxml2=/usr \
-                    --with-magics=/Users/m214003/local/magics-2.14.9 \
 	            CC=gcc CFLAGS="-g -pipe -D_REENTRANT -Wall -W -Wfloat-equal -pedantic -O3 -march=native" \
                     LIBS="-L/opt/local/lib -lopenjpeg"
 #                    --with-libxml2=/usr

src/clipping/clipping.c

@@ -191,7 +191,6 @@ static void compute_lat_circle_z_value(double a[], double b[], double z[]) {
 static unsigned is_inside_gc(double point[], double norm_vec[]) {
 
   double dot;
-  double const tol = 1e-12;
 
   // the product is defined as follows
   // a * b = |a| * |b| * cos(alpha)
@@ -200,7 +199,7 @@ static unsigned is_inside_gc(double point[], double norm_vec[]) {
   dot = dotproduct(point, norm_vec);
 
   // if the point is on the line
-  if (fabs(dot) < tol)
+  if (fabs(dot) < cos(M_PI_2 - angle_tol))
     return 2;
 
   return dot < 0;

src/clipping/geometry.h

@@ -43,6 +43,8 @@
 #include "grid.h"
 #include "utils.h"
 
+extern const double angle_tol;
+
 struct line {
    struct {
       double x;

src/clipping/intersection.c

@@ -41,6 +41,8 @@
 #include "utils.h"
 #include "geometry.h"
 
+// angle tolerance
+const double angle_tol = 1e-10;
 static double const tol = 1.0e-10;
 
 static void crossproduct (double a[], double b[], double cross[]) {

src/remaplib.c

@@ -982,9 +982,9 @@ void remap_vars_init(int map_type, long src_grid_size, long tgt_grid_size, remap
     {
       if      ( map_type == MAP_TYPE_CONSERV   ) rv->sort_add = TRUE;
       else if ( map_type == MAP_TYPE_CONSPHERE ) rv->sort_add = TRUE;
-      else if ( map_type == MAP_TYPE_BILINEAR  ) rv->sort_add = FALSE;
-      else if ( map_type == MAP_TYPE_BICUBIC   ) rv->sort_add = FALSE;
-      else if ( map_type == MAP_TYPE_DISTWGT   ) rv->sort_add = FALSE;
+      else if ( map_type == MAP_TYPE_BILINEAR  ) rv->sort_add = TRUE;
+      else if ( map_type == MAP_TYPE_BICUBIC   ) rv->sort_add = TRUE;
+      else if ( map_type == MAP_TYPE_DISTWGT   ) rv->sort_add = TRUE;
       else cdoAbort("Unknown mapping method!");
     }
   else